Index: LMDZ4/trunk/libf/phylmd/calcul_divers.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/calcul_divers.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/calcul_divers.h	(revision 996)
@@ -12,6 +12,6 @@
 c surface terre
        DO i=1, klon
-         IF(pctsrf_new(i,is_ter).GT.0.) THEN
-            paire_ter(i)=airephy(i)*pctsrf_new(i,is_ter)
+         IF(pctsrf(i,is_ter).GT.0.) THEN
+            paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
          ENDIF 
        ENDDO
Index: LMDZ4/trunk/libf/phylmd/change_srf_frac_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 996)
+++ LMDZ4/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 996)
@@ -0,0 +1,157 @@
+!
+! $Header$
+!
+MODULE change_srf_frac_mod
+
+  IMPLICIT NONE
+
+CONTAINS
+! 
+! Change Surface Fractions
+!
+  SUBROUTINE change_srf_frac(itime, dtime, jour, &
+       pctsrf, alb1, alb2, tsurf, u10m, v10m, pbl_tke)
+!
+! This subroutine is called from physiq.F at each timestep. 
+! 1- For each type of ocean (force, slab, couple) receive new fractions only if
+!    it's time to modify (is_modified=true). Otherwise nothing is done (is_modified=false).   
+! If received new fraction :
+! 2- Tests and ajustements are done on the fractions 
+! 3- Initialize variables where a new fraction(new or melted ice) has appered, 
+!
+
+    USE dimphy 
+    USE surface_data, ONLY : type_ocean
+    USE limit_read_mod
+    USE pbl_surface_mod, ONLY : pbl_surface_newfrac
+    USE cpl_mod, ONLY : cpl_receive_frac
+    USE ocean_slab_mod, ONLY : ocean_slab_frac
+
+    INCLUDE "iniprint.h"
+    INCLUDE "indicesol.h"
+    INCLUDE "YOMCST.h"
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                     :: itime   ! current time step
+    INTEGER, INTENT(IN)                     :: jour    ! day of the year
+    REAL,    INTENT(IN)                     :: dtime   ! length of time step (s)
+  
+! In-Output arguments
+!****************************************************************************************
+   
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf ! sub-surface fraction
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1   ! albedo first interval in SW spektrum
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb2   ! albedo second interval in SW spektrum
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: u10m
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: v10m
+    REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: pbl_tke
+
+! Loccal variables
+!****************************************************************************************
+    INTEGER                        :: i, nsrf
+    LOGICAL                        :: is_modified   ! true if pctsrf is modified at this time step
+    LOGICAL                        :: test_sum=.FALSE.
+    LOGICAL, DIMENSION(klon,nbsrf) :: new_surf
+    REAL, DIMENSION(klon,nbsrf)    :: pctsrf_old    ! fraction from previous time-step
+    REAL                           :: tmpsum
+
+    pctsrf_old(:,:) = pctsrf(:,:)
+!****************************************************************************************
+! 1) 
+! For each type of ocean (force, slab, couple) receive new fractions only if it's time  
+! to modify (is_modified=true). Otherwise nothing is done (is_modified=false).   
+!****************************************************************************************
+    SELECT CASE (type_ocean)
+    CASE ('force')
+       ! Read fraction from limit.nc
+       CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
+    CASE ('slab')
+       ! Get fraction from slab module
+       CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
+    CASE ('couple')
+       ! Get fraction from the coupler
+       CALL cpl_receive_frac(itime, dtime, pctsrf, is_modified)
+    END SELECT
+
+    IF (is_modified) THEN
+!****************************************************************************************
+! 2) 
+! Tests and ajustements on the new fractions :
+! - Put to zero fractions that are too small
+! - Test total fraction sum is one for each grid point
+!
+!****************************************************************************************
+  
+! Test and exit if a fraction is negative
+       IF (MINVAL(pctsrf(:,:)) < 0.) THEN
+          WRITE(lunout,*)'Warning! One or several fractions are negative, itime=',itime
+          WRITE(lunout,*)'at point = ',MINLOC(pctsrf(:,:))
+          WRITE(lunout,*)'value = ',MINVAL(pctsrf(:,:)) 
+          CALL abort_gcm('change_srf_frac','Negative fraction',1)
+       END IF
+
+! Optional test on the incoming fraction 
+       IF (test_sum) THEN
+          DO i= 1, klon
+             tmpsum = SUM(pctsrf(i,:))
+             IF (ABS(1. - tmpsum) > 0.05) CALL abort_gcm('change_srf_frac','Total fraction not equal 1.',1)
+          END DO
+       END IF
+
+! Test for too small fractions of the sum land+landice and ocean+sea-ice
+       WHERE ((pctsrf(:,is_ter) + pctsrf(:,is_lic)) < 2*EPSFRA)
+          pctsrf(:,is_ter) = 0.
+          pctsrf(:,is_lic) = 0.
+       END WHERE
+
+       WHERE ((pctsrf(:,is_oce) + pctsrf(:,is_sic)) < 2*EPSFRA)
+          pctsrf(:,is_oce) = 0.
+          pctsrf(:,is_sic) = 0.
+       END WHERE
+
+! Normalize to force total fraction to be equal one
+       DO i= 1, klon
+          tmpsum = SUM(pctsrf(i,:))
+          DO nsrf = 1, nbsrf
+             pctsrf(i,nsrf) = pctsrf(i,nsrf) / tmpsum
+          END DO
+       END DO
+
+! Test for too small fractions at each sub-surface
+       WHERE (pctsrf(:,is_ter) < EPSFRA)
+          pctsrf(:,is_lic) = pctsrf(:,is_ter) + pctsrf(:,is_lic)
+          pctsrf(:,is_ter) = 0.
+       END WHERE
+
+       WHERE (pctsrf(:,is_lic) < EPSFRA)
+          pctsrf(:,is_ter) = pctsrf(:,is_ter) + pctsrf(:,is_lic)
+          pctsrf(:,is_lic) = 0.
+       END WHERE
+
+       WHERE (pctsrf(:,is_oce) < EPSFRA)
+          pctsrf(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
+          pctsrf(:,is_oce) = 0.
+       END WHERE
+
+       WHERE (pctsrf(:,is_sic) < EPSFRA)
+          pctsrf(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
+          pctsrf(:,is_sic) = 0.
+       END WHERE
+
+!****************************************************************************************
+! 3)
+! Initialize variables where a new fraction has appered, 
+! i.e. where new sea ice has been formed
+! or where ice free ocean has appread in a grid cell
+! 
+!****************************************************************************************
+       CALL pbl_surface_newfrac(itime, pctsrf, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, pbl_tke)
+
+    END IF ! is_modified
+
+  END SUBROUTINE change_srf_frac
+
+
+END MODULE change_srf_frac_mod
Index: LMDZ4/trunk/libf/phylmd/clesphys.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/clesphys.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/clesphys.h	(revision 996)
@@ -47,5 +47,4 @@
        REAL freq_ISCCP, ecrit_ISCCP
        INTEGER :: ip_ebil_phy
-       LOGICAL ok_slab_sicOBS
 
        COMMON/clesphys/cycle_diurne, soil_model, new_oliq,              &
@@ -61,5 +60,5 @@
      &     , ecrit_mth, ecrit_tra, ecrit_reg                            &
      &     , freq_ISCCP, ecrit_ISCCP, ip_ebil_phy                       &
-     &     , ok_slab_sicOBS, ok_lic_melt, cvl_corr                      &
+     &     , ok_lic_melt, cvl_corr                                      &
      &     , qsol0
      
Index: LMDZ4/trunk/libf/phylmd/climb_hq_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/climb_hq_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/climb_hq_mod.F90	(revision 996)
@@ -1,2 +1,5 @@
+!
+! $Header$
+!
 MODULE climb_hq_mod
 !
@@ -274,5 +277,5 @@
     REAL, DIMENSION(klon,klev)               :: h_new, q_new
     REAL, DIMENSION(klon)                    :: psref         
-    INTEGER                                  :: k, i
+    INTEGER                                  :: k, i, ierr
 
 !****************************************************************************************
@@ -350,8 +353,10 @@
 !
 !****************************************************************************************
-    DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H)    
-    DEALLOCATE(gamaq, gamah)
-    DEALLOCATE(Kcoefhq)
-
+    DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)    
+    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
+    DEALLOCATE(gamaq, gamah,stat=ierr)
+    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
+    DEALLOCATE(Kcoefhq,stat=ierr)
+    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
     
   END SUBROUTINE climb_hq_up
Index: LMDZ4/trunk/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/conf_phys.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/conf_phys.F90	(revision 996)
@@ -5,5 +5,5 @@
 !
 
-  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, ok_hf, &
+  subroutine conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, &
  &                     solarlong0,seuil_inversion, &
  &                     fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
@@ -16,5 +16,5 @@
 
    use IOIPSL
-!!!!   USE surface_data,     ONLY : ocean, ok_veget
+   USE surface_data
 
    implicit none
@@ -51,5 +51,5 @@
 ! Sortie:
   character (len = 6)  :: ocean
-  logical              :: ok_veget, ok_newmicro
+  logical              :: ok_newmicro
   integer              :: iflag_radia
   logical              :: ok_journe, ok_mensuel, ok_instan, ok_hf
@@ -116,5 +116,4 @@
   REAL,SAVE :: ecrit_hf_omp, ecrit_day_omp, ecrit_mth_omp, ecrit_reg_omp
   REAL,SAVE :: ecrit_tra_omp
-  LOGICAL, SAVE :: ok_slab_sicOBS_omp
   REAL,SAVE :: cvl_corr_omp
   LOGICAL,SAVE :: ok_lic_melt_omp
@@ -997,14 +996,4 @@
   call getin('ecrit_reg',ecrit_reg_omp)
 !
-!
-!
-!Config Key  = ok_slab_sicOBS
-!Config Desc =
-!Config Def  = .true.
-!Config Help = Pour faire tourner le slab avec fraction 
-!              de glace de mer Observee
-!
-  ok_slab_sicOBS_omp = .true.
-  call getin('ok_slab_sicOBS', ok_slab_sicOBS_omp)
 !
 !
@@ -1198,5 +1187,4 @@
     ecrit_tra = ecrit_tra_omp
     ecrit_reg = ecrit_reg_omp
-    ok_slab_sicOBS = ok_slab_sicOBS_omp
     cvl_corr = cvl_corr_omp
     ok_lic_melt = ok_lic_melt_omp
@@ -1213,7 +1201,27 @@
 !$OMP MASTER
 
+! Attribution of new parmeters according to parameters in .def
+    IF (ocean=='couple' .OR. ocean=='opa8') THEN
+       type_ocean='couple'
+       version_ocean='opa8'
+    ELSE IF (ocean=='nemo') THEN
+       type_ocean='couple'
+       version_ocean='nemo'
+    ELSE IF (ocean=='force') THEN
+       type_ocean='force'
+       version_ocean='xxxxxx'
+    ELSE IF (ocean=='slab') THEN
+       type_ocean='slab'
+       version_ocean='sicOBS'
+    ELSE
+       WRITE(numout,*)' ERROR ocean not valid : ocean= ', ocean
+       CALL abort_gcm('conf_phys','ocean not valid',1)
+    END IF
+
   write(numout,*)' ##############################################'
   write(numout,*)' Configuration des parametres de la physique: '
   write(numout,*)' Config ocean = ', ocean
+  write(numout,*)' Type ocean = ', type_ocean
+  write(numout,*)' Version ocean = ', version_ocean
   write(numout,*)' Config veget = ', ok_veget
   write(numout,*)' Sortie journaliere = ', ok_journe
Index: LMDZ4/trunk/libf/phylmd/cpl_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/cpl_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/cpl_mod.F90	(revision 996)
@@ -32,6 +32,6 @@
   PRIVATE
 
-  ! All subroutine are public except cpl_send_all and cpl_receive_all
-  PUBLIC :: cpl_init, cpl_receive_ocean_fields, cpl_receive_seaice_fields, &
+  ! All subroutine are public except cpl_send_all
+  PUBLIC :: cpl_init, cpl_receive_frac, cpl_receive_ocean_fields, cpl_receive_seaice_fields, &
        cpl_send_ocean_fields, cpl_send_seaice_fields, cpl_send_land_fields, &
        cpl_send_landice_fields, gath2cpl
@@ -68,7 +68,4 @@
   !$OMP THREADPRIVATE(read_alb_sic)
   
-! fraction for different surface, saved during whole coupling period
-  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: pctsrf_sav   
-  !$OMP THREADPRIVATE(pctsrf_sav)
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: unity
   !$OMP THREADPRIVATE(unity)
@@ -88,6 +85,5 @@
   !$OMP THREADPRIVATE(cpl_windsp2D)
   
-! variable for OPENMP parallelisation
-
+! variables for OPENMP parallelisation
   INTEGER,ALLOCATABLE,DIMENSION(:),SAVE :: knon_omp
   REAL,ALLOCATABLE,DIMENSION(:,:),SAVE ::  buffer_omp
@@ -146,6 +142,4 @@
     ALLOCATE(unity(klon), stat = error)
     sum_error = sum_error + error
-    ALLOCATE(pctsrf_sav(klon,nbsrf), stat = error)
-    sum_error = sum_error + error
     ALLOCATE(cpl_sols(klon,2), stat = error) 
     sum_error = sum_error + error
@@ -196,10 +190,9 @@
        unity(ig) = ig
     ENDDO
-    pctsrf_sav = 0.
-
-    cpl_sols = 0.   ; cpl_nsol = 0.  ; cpl_rain = 0.   ; cpl_snow = 0.
-    cpl_evap = 0.   ; cpl_tsol = 0.  ; cpl_fder = 0.   ; cpl_albe = 0.
-    cpl_taux = 0.   ; cpl_tauy = 0.  ; cpl_rriv2D = 0. ; cpl_rcoa2D = 0.
-    cpl_rlic2D = 0. ; cpl_windsp = 0.
+
+!    cpl_sols = 0.   ; cpl_nsol = 0.  ; cpl_rain = 0.   ; cpl_snow = 0.
+!    cpl_evap = 0.   ; cpl_tsol = 0.  ; cpl_fder = 0.   ; cpl_albe = 0.
+!    cpl_taux = 0.   ; cpl_tauy = 0.  ; cpl_rriv2D = 0. ; cpl_rcoa2D = 0.
+!    cpl_rlic2D = 0. ; cpl_windsp = 0.
 
 !*************************************************************************************
@@ -262,6 +255,6 @@
 
 !$OMP MASTER
-  ALLOCATE(knon_omp(0:omp_size-1)) 
-  ALLOCATE(buffer_omp(klon_mpi,0:omp_size-1))       
+    ALLOCATE(knon_omp(0:omp_size-1)) 
+    ALLOCATE(buffer_omp(klon_mpi,0:omp_size-1))       
 !$OMP END MASTER
 !$OMP BARRIER
@@ -272,8 +265,11 @@
 !*************************************************************************************
 !
-
-  SUBROUTINE cpl_receive_all(itime, dtime, pctsrf)
-! This subroutine reads from coupler for both ocean and seaice
+ 
+  SUBROUTINE cpl_receive_frac(itime, dtime, pctsrf, is_modified)
+! This subroutine receives from coupler for both ocean and seaice
 ! 4 fields : read_sst, read_sic, read_sit and read_alb_sic. 
+! The new sea-ice-land-landice fraction is returned. The others fields 
+! are stored in this module.
+    USE surface_data
 
     INCLUDE "indicesol.h"
@@ -283,20 +279,21 @@
     INCLUDE "dimensions.h"
 
-! Input arguments
+! Arguments
 !************************************************************************************
-    INTEGER, INTENT(IN)                     :: itime
-    REAL, INTENT(IN)                        :: dtime
-    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
+    INTEGER, INTENT(IN)                        :: itime
+    REAL, INTENT(IN)                           :: dtime
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf
+    LOGICAL, INTENT(OUT)                       :: is_modified
 
 ! Local variables
 !************************************************************************************
-    INTEGER                                 :: j, ig, il_time_secs
+    INTEGER                                 :: j, i, time_sec
     INTEGER                                 :: itau_w
     INTEGER, DIMENSION(iim*(jjm+1))         :: ndexcs
-    CHARACTER(len = 20)                     :: modname = 'cpl_receive_all'
+    CHARACTER(len = 20)                     :: modname = 'cpl_receive_frac'
     CHARACTER(len = 80)                     :: abort_message
     REAL, DIMENSION(klon)                   :: read_sic1D
-    REAL, DIMENSION(iim,jj_nb,jpfldo2a)  :: tab_read_flds
-    REAL, DIMENSION(iim,jj_nb)           :: read_sic
+    REAL, DIMENSION(iim,jj_nb,jpfldo2a)     :: tab_read_flds
+    REAL, DIMENSION(klon,nbsrf)             :: pctsrf_old
 
 !*************************************************************************************
@@ -306,117 +303,72 @@
 !*************************************************************************************
 
+    is_modified=.FALSE.
+
+! Check if right moment to recevie from coupler
+    IF (MOD(itime, nexca) == 1) THEN
+       is_modified=.TRUE.
+ 
+       time_sec=(itime-1)*dtime
 #ifdef CPP_COUPLE
-    il_time_secs=(itime-1)*dtime
 !$OMP MASTER
-    CALL fromcpl(il_time_secs, tab_read_flds)
+       CALL fromcpl(time_sec, tab_read_flds)
 !$OMP END MASTER
 #endif
     
 ! NetCDF output of received fields
-    IF (is_sequential) THEN
-       ndexcs(:) = 0
-       itau_w = itau_phy + itime
-       DO ig = 1, jpfldo2a
-          CALL histwrite(nidcs,cl_read(ig),itau_w,tab_read_flds(:,:,ig),iim*(jjm+1),ndexcs)
-       END DO
-    ENDIF
-
+       IF (is_sequential) THEN
+          ndexcs(:) = 0
+          itau_w = itau_phy + itime
+          DO i = 1, jpfldo2a
+             CALL histwrite(nidcs,cl_read(i),itau_w,tab_read_flds(:,:,i),iim*(jjm+1),ndexcs)
+          END DO
+       ENDIF
+
+! Save each field in a 2D array. 
 !$OMP MASTER
-
-! Save each field in a 2D array. 
-
-    IF (OPA_version=='OPA9') THEN
-      read_sst(:,:)     = tab_read_flds(:,:,1)  ! Sea surface temperature
-      read_sic(:,:)     = tab_read_flds(:,:,2)  ! Sea ice concentration
-      read_sit(:,:)     = tab_read_flds(:,:,3)  ! Sea ice temperature
-      read_alb_sic(:,:) = tab_read_flds(:,:,4)  ! Albedo at sea ice
-    ELSE IF (OPA_version=='OPA8') THEN
-      read_sst(:,:)     = tab_read_flds(:,:,1)  ! Sea surface temperature
-      read_sic(:,:)     = tab_read_flds(:,:,2)  ! Sea ice concentration
-      read_alb_sic(:,:) = tab_read_flds(:,:,3)  ! Albedo at sea ice
-      read_sit(:,:)     = tab_read_flds(:,:,4)  ! Sea ice temperature
-    ELSE
-      STOP 'Bad OPA version for coupled model'
-    ENDIF
-
-!*************************************************************************************
-! Temperature and albedo are weighted with the fraction of sea-ice(read-sic)
-!
-!*************************************************************************************
-    DO j = 1, jj_nb
-       DO ig = 1, iim
-          IF (ABS(1. - read_sic(ig,j)) < 0.00001) THEN
-             read_sst(ig,j) = RTT - 1.8
-             read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
-             read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
-          ELSE IF (ABS(read_sic(ig,j)) < 0.00001) THEN
-             read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
-             read_sit(ig,j) = read_sst(ig,j)
-             read_alb_sic(ig,j) =  0.6
-          ELSE
-             read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
-             read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
-             read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
+       IF (version_ocean=='nemo') THEN
+          read_sst(:,:)     = tab_read_flds(:,:,1)  ! Sea surface temperature
+          read_sic(:,:)     = tab_read_flds(:,:,2)  ! Sea ice concentration
+          read_sit(:,:)     = tab_read_flds(:,:,3)  ! Sea ice temperature
+          read_alb_sic(:,:) = tab_read_flds(:,:,4)  ! Albedo at sea ice
+       ELSE IF (version_ocean=='opa8') THEN
+          read_sst(:,:)     = tab_read_flds(:,:,1)  ! Sea surface temperature (multiplicated by fraction)
+          read_sic(:,:)     = tab_read_flds(:,:,2)  ! Sea ice concentration 
+          read_alb_sic(:,:) = tab_read_flds(:,:,3)  ! Albedo at sea ice (multiplicated by fraction)
+          read_sit(:,:)     = tab_read_flds(:,:,4)  ! Sea ice temperature (multiplicated by fraction)
+       END IF
+!$OMP END MASTER
+
+!*************************************************************************************
+!  Transform seaice fraction (read_sic : ocean-seaice mask) into global 
+!  fraction (pctsrf : ocean-seaice-land-landice mask)
+!
+!*************************************************************************************
+       CALL cpl2gath(read_sic, read_sic1D, klon, unity)
+
+       pctsrf_old(:,:) = pctsrf(:,:)
+       DO i = 1, klon
+          ! treatment only of points with ocean and/or seaice
+          IF (pctsrf_old(i,is_oce) + pctsrf_old(i,is_sic) > 0.) THEN
+             pctsrf(i,is_sic) = (pctsrf_old(i,is_oce) + pctsrf_old(i,is_sic)) &
+                  * read_sic1D(i)
+             pctsrf(i,is_oce) = (pctsrf_old(i,is_oce) + pctsrf_old(i,is_sic)) &
+                  - pctsrf(i,is_sic)
           ENDIF
        ENDDO
-    ENDDO
-!$OMP END MASTER
-
-!*************************************************************************************
-!  Transform seaice fraction, read_sic into pctsrf_sav
-!
-!*************************************************************************************
-    CALL cpl2gath(read_sic, read_sic1D, klon, unity)
-
-    DO ig = 1, klon
-       ! treatment only of ocean and/or seaice points
-       IF (pctsrf(ig,is_oce) > epsfra .OR. &
-            pctsrf(ig,is_sic) > epsfra) THEN
-          pctsrf_sav(ig,is_sic) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
-               * read_sic1D(ig)
-          pctsrf_sav(ig,is_oce) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
-               - pctsrf_sav(ig,is_sic)
-       ENDIF
-    ENDDO
-
-!*************************************************************************************
-! To avoid round up problems
-!
-!*************************************************************************************
-    WHERE (ABS(pctsrf_sav(:,is_sic)) .LE. 2.*EPSILON(pctsrf_sav(1,is_sic)))
-       pctsrf_sav(:,is_sic) = 0.
-       pctsrf_sav(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
-    ENDWHERE
-    WHERE (ABS(pctsrf_sav(:,is_oce)) .LE. 2.*EPSILON(pctsrf_sav(1,is_oce)))
-       pctsrf_sav(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
-       pctsrf_sav(:,is_oce) = 0.
-    ENDWHERE
-    IF (MINVAL(pctsrf_sav(:,is_oce)) < 0.) THEN
-       WRITE(*,*)'Pb fraction ocean inferieure a 0'
-       WRITE(*,*)'au point ',MINLOC(pctsrf_sav(:,is_oce))
-       WRITE(*,*)'valeur = ',MINVAL(pctsrf_sav(:,is_oce)) 
-       abort_message = 'voir ci-dessus'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-    IF (MINVAL(pctsrf_sav(:,is_sic)) < 0.) THEN
-       WRITE(*,*)'Pb fraction glace inferieure a 0'
-       WRITE(*,*)'au point ',MINLOC(pctsrf_sav(:,is_sic))
-       WRITE(*,*)'valeur = ',MINVAL(pctsrf_sav(:,is_sic)) 
-       abort_message = 'voir ci-dessus'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-        
-  END SUBROUTINE cpl_receive_all
-!
-!*************************************************************************************
-! 
-  SUBROUTINE cpl_receive_ocean_fields(itime, dtime, knon, knindex, pctsrf, &
-       tsurf_new, pctsrf_oce)
-!
-! This routine reads, if first time step in a coupling period, all fields reveived from
-! coupler for all types of surfaces. It returns the fields for the ocean surface which
-! are the sea surface temperature and the fraction of ocean. 
-! The fields are transformed into 1D arrays with valid points : 
-! tsurf_new(1:knon), pctsrf(1:klon).
+
+    END IF ! if time to receive
+
+  END SUBROUTINE cpl_receive_frac
+
+!
+!*************************************************************************************
+!
+
+  SUBROUTINE cpl_receive_ocean_fields(knon, knindex, tsurf_new)
+!
+! This routine returns the field for the ocean that has been read from the coupler
+! (done earlier with cpl_receive_frac). The field is the temperature.
+! The temperature is transformed into 1D array with valid points from index 1 to knon.
 !
     INCLUDE "indicesol.h"
@@ -424,23 +376,15 @@
 ! Input arguments
 !*************************************************************************************
-    INTEGER, INTENT(IN)                     :: itime
-    REAL, INTENT(IN)                        :: dtime
     INTEGER, INTENT(IN)                     :: knon
     INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
-    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
 
 ! Output arguments
 !*************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)      :: tsurf_new
-    REAL, DIMENSION(klon), INTENT(OUT)      :: pctsrf_oce
-
-!*************************************************************************************
-! If first time step in a coupling period receive all fields for all types
-! of surfaces from coupler : read_sst, read_sit, read_alb_sic and pctsrf_sav. 
-! 
-!*************************************************************************************
-
-    IF (MOD(itime, nexca) == 1) CALL cpl_receive_all(itime, dtime, pctsrf)
-
+
+! Local variables
+!*************************************************************************************
+    INTEGER               :: i
+    REAL, DIMENSION(klon) :: sic_new
 
 !*************************************************************************************
@@ -449,21 +393,29 @@
 !*************************************************************************************
     CALL cpl2gath(read_sst, tsurf_new, knon, knindex)
-    pctsrf_oce(:) = pctsrf_sav(:,is_oce) 
-
+    CALL cpl2gath(read_sic, sic_new, knon, knindex)
+
+!*************************************************************************************
+! The fields received from the coupler have to be weighted with the fraction of ocean 
+! in relation to the total sea-ice+ocean
+!
+!*************************************************************************************
+    DO i=1, knon
+       tsurf_new(i) = tsurf_new(i)/(1. - sic_new(i))
+    END DO
 
   END SUBROUTINE cpl_receive_ocean_fields
-!
-!*************************************************************************************
-!
+
+!
+!*************************************************************************************
+!
+
   SUBROUTINE cpl_receive_seaice_fields(knon, knindex, &
-       tsurf_new, alb_new, pctsrf_sic)
+       tsurf_new, alb_new)
 !
 ! This routine returns the fields for the seaice that have been read from the coupler
-! (done earlier with cpl_receive_ocean_fields). These fields are the temperature and 
+! (done earlier with cpl_receive_frac). These fields are the temperature and 
 ! albedo at sea ice surface and fraction of sea ice.
-! The fields are transformed into  1D arrays with valid points : 
-! tsurf_new(1:knon), alb_new(1:knon), pctsrf(1:klon).
-!
-    INCLUDE "indicesol.h"
+! The fields are transformed into 1D arrays with valid points from index 1 to knon. 
+!
 
 ! Input arguments
@@ -476,6 +428,9 @@
     REAL, DIMENSION(klon), INTENT(OUT)      :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)      :: alb_new
-    REAL, DIMENSION(klon), INTENT(OUT)      :: pctsrf_sic
-
+
+! Local variables
+!*************************************************************************************
+    INTEGER               :: i
+    REAL, DIMENSION(klon) :: sic_new
 
 !*************************************************************************************
@@ -485,5 +440,15 @@
     CALL cpl2gath(read_sit, tsurf_new, knon, knindex)
     CALL cpl2gath(read_alb_sic, alb_new, knon, knindex)
-    pctsrf_sic(:) = pctsrf_sav(:,is_sic) 
+    CALL cpl2gath(read_sic, sic_new, knon, knindex)
+
+!*************************************************************************************
+! The fields received from the coupler have to be weighted with the sea-ice 
+! concentration (in relation to the total sea-ice + ocean).
+!
+!*************************************************************************************
+    DO i= 1, knon
+       tsurf_new(i) = tsurf_new(i) / sic_new(i)
+       alb_new(i)   = alb_new(i)   / sic_new(i)
+    END DO
 
   END SUBROUTINE cpl_receive_seaice_fields
@@ -535,15 +500,15 @@
 !*************************************************************************************
     IF (MOD(itime, nexca) == 1) THEN
-       cpl_sols(:,cpl_index) = 0.0
-       cpl_nsol(:,cpl_index) = 0.0
-       cpl_rain(:,cpl_index) = 0.0
-       cpl_snow(:,cpl_index) = 0.0
-       cpl_evap(:,cpl_index) = 0.0
-       cpl_tsol(:,cpl_index) = 0.0
-       cpl_fder(:,cpl_index) = 0.0
-       cpl_albe(:,cpl_index) = 0.0
-       cpl_taux(:,cpl_index) = 0.0
-       cpl_tauy(:,cpl_index) = 0.0
-       cpl_windsp(:,cpl_index) = 0.0
+       cpl_sols(1:knon,cpl_index) = 0.0
+       cpl_nsol(1:knon,cpl_index) = 0.0
+       cpl_rain(1:knon,cpl_index) = 0.0
+       cpl_snow(1:knon,cpl_index) = 0.0
+       cpl_evap(1:knon,cpl_index) = 0.0
+       cpl_tsol(1:knon,cpl_index) = 0.0
+       cpl_fder(1:knon,cpl_index) = 0.0
+       cpl_albe(1:knon,cpl_index) = 0.0
+       cpl_taux(1:knon,cpl_index) = 0.0
+       cpl_tauy(1:knon,cpl_index) = 0.0
+       cpl_windsp(1:knon,cpl_index) = 0.0
     ENDIF
        
@@ -709,14 +674,14 @@
 !*************************************************************************************
     IF (MOD(itime, nexca) == 1) THEN
-       cpl_sols(:,cpl_index) = 0.0
-       cpl_nsol(:,cpl_index) = 0.0
-       cpl_rain(:,cpl_index) = 0.0
-       cpl_snow(:,cpl_index) = 0.0
-       cpl_evap(:,cpl_index) = 0.0
-       cpl_tsol(:,cpl_index) = 0.0
-       cpl_fder(:,cpl_index) = 0.0
-       cpl_albe(:,cpl_index) = 0.0
-       cpl_taux(:,cpl_index) = 0.0
-       cpl_tauy(:,cpl_index) = 0.0
+       cpl_sols(1:knon,cpl_index) = 0.0
+       cpl_nsol(1:knon,cpl_index) = 0.0
+       cpl_rain(1:knon,cpl_index) = 0.0
+       cpl_snow(1:knon,cpl_index) = 0.0
+       cpl_evap(1:knon,cpl_index) = 0.0
+       cpl_tsol(1:knon,cpl_index) = 0.0
+       cpl_fder(1:knon,cpl_index) = 0.0
+       cpl_albe(1:knon,cpl_index) = 0.0
+       cpl_taux(1:knon,cpl_index) = 0.0
+       cpl_tauy(1:knon,cpl_index) = 0.0
     ENDIF
        
@@ -893,5 +858,5 @@
 !
     INCLUDE "dimensions.h"
-    
+
 ! Input varibales
 !*************************************************************************************
@@ -944,4 +909,5 @@
 ! all calculations at the different surfaces have to be done before. 
 !    
+    USE surface_data
 ! Some includes
 !*************************************************************************************
@@ -962,5 +928,5 @@
     INTEGER                                              :: error, sum_error, j
     INTEGER                                              :: itau_w
-    INTEGER                                              :: il_time_secs
+    INTEGER                                              :: time_sec
     INTEGER, DIMENSION(iim*(jjm+1))                      :: ndexct
     REAL                                                 :: Up, Down
@@ -994,32 +960,28 @@
 !
 !*************************************************************************************
-!! AC >>
-
 !$OMP MASTER
-    IF (OPA_version=='OPA9') THEN
-      tab_flds(:,:,7)  = cpl_windsp2D(:,:)
-      tab_flds(:,:,14) = cpl_sols2D(:,:,2)
-      tab_flds(:,:,12) = cpl_sols2D(:,:,1)
-      tab_flds(:,:,15) = cpl_nsol2D(:,:,2)
-      tab_flds(:,:,13) = cpl_nsol2D(:,:,1)
-      tab_flds(:,:,16) = cpl_fder2D(:,:,2)
-      tab_flds(:,:,11) = cpl_evap2D(:,:,2)
-      tab_flds(:,:,18) = cpl_rriv2D(:,:)
-      tab_flds(:,:,19) = cpl_rcoa2D(:,:)
-    ELSE IF (OPA_version=='OPA8') THEN
-      tab_flds(:,:,7)  = cpl_windsp2D(:,:)
-      tab_flds(:,:,8)  = cpl_sols2D(:,:,2)
-      tab_flds(:,:,9)  = cpl_sols2D(:,:,1)
-      tab_flds(:,:,10) = cpl_nsol2D(:,:,2)
-      tab_flds(:,:,11) = cpl_nsol2D(:,:,1)
-      tab_flds(:,:,12) = cpl_fder2D(:,:,2)
-      tab_flds(:,:,13) = cpl_evap2D(:,:,2)
-      tab_flds(:,:,14) = cpl_evap2D(:,:,1)
-      tab_flds(:,:,17) = cpl_rcoa2D(:,:)
-      tab_flds(:,:,18) = cpl_rriv2D(:,:)
-    ELSE
-      STOP 'Bad OPA version for coupled model'
-    ENDIF
-
+    IF (version_ocean=='nemo') THEN
+       tab_flds(:,:,7)  = cpl_windsp2D(:,:)
+       tab_flds(:,:,14) = cpl_sols2D(:,:,2)
+       tab_flds(:,:,12) = cpl_sols2D(:,:,1)
+       tab_flds(:,:,15) = cpl_nsol2D(:,:,2)
+       tab_flds(:,:,13) = cpl_nsol2D(:,:,1)
+       tab_flds(:,:,16) = cpl_fder2D(:,:,2)
+       tab_flds(:,:,11) = cpl_evap2D(:,:,2)
+       tab_flds(:,:,18) = cpl_rriv2D(:,:)
+       tab_flds(:,:,19) = cpl_rcoa2D(:,:)
+    ELSE IF (version_ocean=='opa8') THEN
+       tab_flds(:,:,7)  = cpl_windsp2D(:,:)
+       tab_flds(:,:,8)  = cpl_sols2D(:,:,2)
+       tab_flds(:,:,9)  = cpl_sols2D(:,:,1)
+       tab_flds(:,:,10) = cpl_nsol2D(:,:,2)
+       tab_flds(:,:,11) = cpl_nsol2D(:,:,1)
+       tab_flds(:,:,12) = cpl_fder2D(:,:,2)
+       tab_flds(:,:,13) = cpl_evap2D(:,:,2)
+       tab_flds(:,:,14) = cpl_evap2D(:,:,1)
+       tab_flds(:,:,17) = cpl_rcoa2D(:,:)
+       tab_flds(:,:,18) = cpl_rriv2D(:,:)
+    END IF
+    
 !*************************************************************************************
 ! Transform the fraction of sub-surfaces from 1D to 2D array
@@ -1028,5 +990,4 @@
     pctsrf2D(:,:,:) = 0.
 !$OMP END MASTER
-
     CALL gath2cpl(pctsrf(:,is_oce), pctsrf2D(:,:,is_oce), klon, unity)
     CALL gath2cpl(pctsrf(:,is_sic), pctsrf2D(:,:,is_sic), klon, unity)
@@ -1040,11 +1001,11 @@
     IF (is_omp_root) THEN
 
-      DO j = 1, jj_nb
-         tmp_calv(:,j) = DOT_PRODUCT (cpl_rlic2D(1:iim,j), &
-              pctsrf2D(1:iim,j,is_lic)) / REAL(iim)
-      ENDDO
-    
-    
-      IF (is_parallel) THEN
+       DO j = 1, jj_nb
+          tmp_calv(:,j) = DOT_PRODUCT (cpl_rlic2D(1:iim,j), &
+               pctsrf2D(1:iim,j,is_lic)) / REAL(iim)
+       ENDDO
+       
+    
+       IF (is_parallel) THEN
          IF (.NOT. is_north_pole) THEN
 #ifdef CPP_PARA
@@ -1053,5 +1014,5 @@
 #endif
          ENDIF
-       
+          
          IF (.NOT. is_south_pole) THEN
 #ifdef CPP_PARA
@@ -1060,10 +1021,10 @@
 #endif
          ENDIF
-       
+         
          IF (.NOT. is_north_pole .AND. ii_begin /=1) THEN
             Up=Up+tmp_calv(iim,1)
             tmp_calv(:,1)=Up
          ENDIF
-       
+         
          IF (.NOT. is_south_pole .AND. ii_end /= iim) THEN
             Down=Down+tmp_calv(1,jj_nb)
@@ -1071,13 +1032,10 @@
          ENDIF
       ENDIF
-        
-      IF (OPA_version=='OPA9') THEN
-        tab_flds(:,:,17) = tmp_calv(:,:)
-      ELSE IF (OPA_version=='OPA8') THEN
-        tab_flds(:,:,17) = tmp_calv(:,:)
-      ELSE
-        STOP 'Bad OPA version for coupled model'
-      ENDIF
-
+      
+      IF (version_ocean=='nemo') THEN
+         tab_flds(:,:,17) = tmp_calv(:,:)
+      ELSE IF (version_ocean=='opa8') THEN
+         tab_flds(:,:,19) = tmp_calv(:,:)
+      END IF
 
 !*************************************************************************************
@@ -1089,53 +1047,45 @@
 !
 !*************************************************************************************    
-      ! fraction oce+seaice
-      deno =  pctsrf2D(:,:,is_oce) + pctsrf2D(:,:,is_sic) 
-
-      IF (OPA_version=='OPA9') THEN
-
-        tab_flds(:,:,10) = 0.0
-        tmp_taux(:,:)    = 0.0
-        tmp_tauy(:,:)    = 0.0
-        ! fraction oce+seaice
-        ! For all valid grid cells containing some fraction of ocean or sea-ice
-        WHERE ( deno(:,:) /= 0 )
-           tab_flds(:,:,10) = cpl_snow2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
-                              cpl_snow2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
-
-           tmp_taux = cpl_taux2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
-                      cpl_taux2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
-           tmp_tauy = cpl_tauy2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
-                      cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
-        ENDWHERE
-        tab_flds(:,:,8) = (cpl_evap2D(:,:,1) - ( cpl_rain2D(:,:,1) + cpl_snow2D(:,:,1)))
-        tab_flds(:,:,9) = (cpl_evap2D(:,:,2) - ( cpl_rain2D(:,:,2) + cpl_snow2D(:,:,2)))
-
-      ELSE IF (OPA_version=='OPA8') THEN
-
-        tab_flds(:,:,15) = 0.0
-        tab_flds(:,:,16) = 0.0
-        tmp_taux(:,:)    = 0.0
-        tmp_tauy(:,:)    = 0.0
-        ! For all valid grid cells containing some fraction of ocean or sea-ice
-        WHERE ( deno(:,:) /= 0 )
-           tab_flds(:,:,15) = cpl_rain2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
-                              cpl_rain2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
-           tab_flds(:,:,16) = cpl_snow2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
-                              cpl_snow2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
-       
-           tmp_taux = cpl_taux2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
-                      cpl_taux2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
-           tmp_tauy = cpl_tauy2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
-                      cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
-        ENDWHERE
-
-      ELSE
-        STOP 'Bad OPA version for coupled model'
-      ENDIF
-    
-    ENDIF ! is_omp_root  
-
-
-! AC <<
+       ! fraction oce+seaice
+       deno =  pctsrf2D(:,:,is_oce) + pctsrf2D(:,:,is_sic) 
+
+       IF (version_ocean=='nemo') THEN
+          tab_flds(:,:,10) = 0.0
+          tmp_taux(:,:)    = 0.0
+          tmp_tauy(:,:)    = 0.0
+          ! For all valid grid cells containing some fraction of ocean or sea-ice
+          WHERE ( deno(:,:) /= 0 )
+             tab_flds(:,:,10) = cpl_snow2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
+                  cpl_snow2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
+             
+             tmp_taux = cpl_taux2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
+                  cpl_taux2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
+             tmp_tauy = cpl_tauy2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
+                  cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
+          ENDWHERE
+          tab_flds(:,:,8) = (cpl_evap2D(:,:,1) - ( cpl_rain2D(:,:,1) + cpl_snow2D(:,:,1)))
+          tab_flds(:,:,9) = (cpl_evap2D(:,:,2) - ( cpl_rain2D(:,:,2) + cpl_snow2D(:,:,2)))
+          
+       ELSE IF (version_ocean=='opa8') THEN
+          tab_flds(:,:,15) = 0.0
+          tab_flds(:,:,16) = 0.0
+          tmp_taux(:,:)    = 0.0
+          tmp_tauy(:,:)    = 0.0
+          ! For all valid grid cells containing some fraction of ocean or sea-ice
+          WHERE ( deno(:,:) /= 0 )
+             tab_flds(:,:,15) = cpl_rain2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
+                  cpl_rain2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
+             tab_flds(:,:,16) = cpl_snow2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
+                  cpl_snow2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
+             
+             tmp_taux = cpl_taux2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
+                  cpl_taux2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
+             tmp_tauy = cpl_tauy2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
+                  cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
+          ENDWHERE
+       END IF
+
+    ENDIF ! is_omp_root
+  
 !*************************************************************************************
 ! Transform the wind components from local atmospheric 2D coordinates to geocentric 
@@ -1145,5 +1095,4 @@
 
 ! Transform the longitudes and latitudes on 2D arrays
-    
     CALL gather_omp(rlon,rlon_mpi)
     CALL gather_omp(rlat,rlat_mpi)
@@ -1211,8 +1160,8 @@
 !
 !*************************************************************************************
+    time_sec=(itime-1)*dtime
 #ifdef CPP_COUPLE
-    il_time_secs=(itime-1)*dtime
 !$OMP MASTER
-    CALL intocpl(il_time_secs, lafin, tab_flds(:,:,:))
+    CALL intocpl(time_sec, lafin, tab_flds(:,:,:))
 !$OMP END MASTER
 #endif
@@ -1239,5 +1188,5 @@
 !
   SUBROUTINE cpl2gath(champ_in, champ_out, knon, knindex)
-  USE mod_phys_lmdz_para
+    USE mod_phys_lmdz_para
 ! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
 ! au coupleur.
@@ -1269,6 +1218,4 @@
 !*************************************************************************************
 !
-    
-
 ! Transform from 2 dimensions (iim,jj_nb) to 1 dimension (klon)
 !$OMP MASTER 
@@ -1283,5 +1230,4 @@
        champ_out(i) = temp_omp(ig)
     ENDDO
-  
    
   END SUBROUTINE cpl2gath
Index: LMDZ4/trunk/libf/phylmd/ini_paramLMDZ_phy.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/ini_paramLMDZ_phy.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/ini_paramLMDZ_phy.h	(revision 996)
@@ -31,10 +31,4 @@
      .                "once", zstophy,zout)
 c
-       CALL histdef(nid_ctesGCM, "ok_slab_sicOBS",
-     ."ok_slab_sicOBS= 1: glace de mer observee, =0: gl.mer calculee 
-     . par le slab", "-",
-     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
-     .                "once", zstophy,zout)
-c
        CALL histdef(nid_ctesGCM, "type_run",
      .        "Type run: 1= CLIM ou ENSP, 2= AMIP ou CFMI",
Index: LMDZ4/trunk/libf/phylmd/limit_read_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/limit_read_mod.F90	(revision 996)
+++ LMDZ4/trunk/libf/phylmd/limit_read_mod.F90	(revision 996)
@@ -0,0 +1,331 @@
+!
+! $Header$
+!
+MODULE limit_read_mod
+!
+! This module reads the fichier "limit.nc" containing fields for surface forcing.
+!
+! Module subroutines :
+!  limit_read_frac    : call limit_read_tot and return the fractions
+!  limit_read_rug_alb : return rugosity and albedo, if coupled ocean call limit_read_tot first
+!  limit_read_sst     : return sea ice temperature   
+!  limit_read_tot     : read limit.nc and store the fields in local modules variables
+!
+  IMPLICIT NONE
+
+  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE, PRIVATE :: pctsrf
+  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: rugos
+  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: albedo
+  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: sst
+  LOGICAL :: read_continents=.FALSE.
+
+CONTAINS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!
+!! Public subroutines :
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
+!
+! This subroutine is called from "change_srf_frac" for case of 
+! ocean=force or from ocean_slab_frac for ocean=slab.
+! The fraction for all sub-surfaces at actual time step is returned.
+
+    USE dimphy
+    INCLUDE "indicesol.h"
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN) :: itime   ! time step
+    INTEGER, INTENT(IN) :: jour    ! current day
+    REAL   , INTENT(IN) :: dtime   ! length of time step
+  
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new  ! sub surface fractions
+    LOGICAL, INTENT(OUT)                     :: is_modified ! true if pctsrf is modified at this time step
+
+! End declaration
+!****************************************************************************************
+
+! 1) Read file limit.nc
+    CALL limit_read_tot(itime, dtime, jour, is_modified)
+
+! 2) Return the fraction read in limit_read_tot
+    pctsrf_new(:,:) = pctsrf(:,:)
+    
+  END SUBROUTINE limit_read_frac
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE limit_read_rug_alb(itime, dtime, jour, &
+       knon, knindex, &
+       rugos_out, alb_out)
+!
+! This subroutine is called from surf_land_bucket. 
+! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple"
+! then this routine will call limit_read_tot.
+!
+    USE dimphy
+    USE surface_data
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN) :: itime                     ! numero du pas de temps courant
+    INTEGER, INTENT(IN) :: jour                      ! jour a lire dans l'annee
+    REAL   , INTENT(IN) :: dtime                     ! pas de temps de la physique (en s)
+    INTEGER, INTENT(IN) :: knon                      ! nomber of points on compressed grid
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out
+    REAL, DIMENSION(klon), INTENT(OUT) :: alb_out
+    
+! Local variables
+!****************************************************************************************
+    INTEGER :: i
+    LOGICAL :: is_modified
+!****************************************************************************************
+
+    IF (type_ocean == 'couple') THEN
+       ! limit.nc has not yet been read. Do it now!
+       CALL limit_read_tot(itime, dtime, jour, is_modified)
+    END IF
+
+    DO i=1,knon
+       rugos_out(i) = rugos(knindex(i))
+       alb_out(i)  = albedo(knindex(i))
+    END DO
+
+  END SUBROUTINE limit_read_rug_alb
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE limit_read_sst(knon, knindex, sst_out)
+!
+! This subroutine returns the sea surface temperature already read from limit.nc.
+!
+    USE dimphy, ONLY : klon
+
+    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
+    REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out
+
+    INTEGER :: i
+
+    DO i = 1, knon
+       sst_out(i) = sst(knindex(i))
+    END DO
+
+  END SUBROUTINE limit_read_sst
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!
+!! Private subroutine :
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified)
+!
+! Read everything needed from limit.nc
+!
+! 0) Initialize
+! 1) Open the file limit.nc, if it is time
+! 2) Read fraction, if not type_ocean=couple
+! 3) Read sea surface temperature, if not type_ocean=couple
+! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
+! 5) Close file and distribuate variables to all processus
+
+    USE dimphy
+    USE mod_grid_phy_lmdz
+    USE mod_phys_lmdz_para
+    USE surface_data, ONLY : type_ocean, ok_veget
+    USE netcdf
+
+    IMPLICIT NONE
+    
+    INCLUDE "indicesol.h"
+
+! In- and ouput arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN) :: itime   ! numero du pas de temps courant
+    INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
+    REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
+
+    LOGICAL, INTENT(OUT) :: is_modified  ! true if pctsrf is modified at this time step
+
+! Locals variables with attribute SAVE
+!****************************************************************************************
+! frequence de lecture des conditions limites (en pas de physique) 
+    INTEGER,SAVE                              :: lmt_pas   
+    LOGICAL, SAVE                             :: first_call=.TRUE.
+    
+! Locals variables
+!****************************************************************************************
+    INTEGER                                   :: nid, nvarid
+    INTEGER                                   :: ii, ierr
+    INTEGER, DIMENSION(2)                     :: start, epais
+    REAL, DIMENSION(klon_glo,nbsrf)           :: pct_glo  ! fraction at global grid
+    REAL, DIMENSION(klon_glo)                 :: sst_glo  ! sea-surface temperature at global grid
+    REAL, DIMENSION(klon_glo)                 :: rug_glo  ! rugosity at global grid
+    REAL, DIMENSION(klon_glo)                 :: alb_glo  ! albedo at global grid
+    CHARACTER(len=20)                         :: modname='limit_read_mod'     
+
+! End declaration
+!****************************************************************************************
+
+!****************************************************************************************
+! 0) Initialization
+!
+!****************************************************************************************
+    IF (first_call) THEN
+       ! calculate number of time steps for one day
+       lmt_pas = NINT(86400./dtime * 1.0)
+       
+       ! Allocate module save variables
+       IF ( type_ocean /= 'couple' ) THEN
+          ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
+          IF (ierr /= 0) CALL abort_gcm(modname, 'PB in allocating pctsrf and sst',1)
+       END IF
+
+       IF ( .NOT. ok_veget ) THEN
+          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
+          IF (ierr /= 0) CALL abort_gcm(modname, 'PB in allocating rugos and albedo',1)
+       END IF
+
+       first_call=.FALSE.
+    ENDIF
+  
+!****************************************************************************************
+! 1) Open the file limit.nc if it is the right moment to read, once a day.
+!    The file is read only by the master thread of the master mpi process(is_mpi_root)
+!
+!****************************************************************************************
+
+    is_modified = .FALSE.
+    IF (MOD(itime-1, lmt_pas) == 0) THEN   ! time to read
+       is_modified = .TRUE.
+!$OMP MASTER  ! Only master thread
+       IF (is_mpi_root) THEN ! Only master processus
+
+          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
+          IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,&
+               'Pb d''ouverture du fichier de conditions aux limites',1)
+          
+          ! La tranche de donnees a lire:
+          start(1) = 1
+          start(2) = jour
+          epais(1) = klon_glo
+          epais(2) = 1
+
+
+!****************************************************************************************
+! 2) Read fraction if not type_ocean=couple
+!
+!****************************************************************************************
+
+          IF ( type_ocean /= 'couple') THEN
+!
+! Ocean fraction
+             ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname, 'Le champ <FOCE> est absent',1)
+             
+             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <FOCE>' ,1)
+!
+! Sea-ice fraction
+             ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <FSIC> est absent',1)
+
+             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <FSIC>' ,1)
+
+
+! Read land and continentals fraction only if asked for
+             IF (read_continents .OR. itime == 1) THEN
+!
+! Land fraction
+                ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
+                IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <FTER> est absent',1)
+                
+                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
+                IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <FTER>',1)
+!
+! Continentale ice fraction
+                ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
+                IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <FLIC> est absent',1)
+
+                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
+                IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <FLIC>',1)
+             END IF
+
+          END IF ! type_ocean /= couple
+
+!****************************************************************************************
+! 3) Read sea-surface temperature, if not coupled ocean
+!
+!****************************************************************************************
+          IF ( type_ocean /= 'couple') THEN
+
+             ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <SST> est absent',1)
+
+             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <SST>',1)
+          
+          END IF
+
+!****************************************************************************************
+! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
+!
+!****************************************************************************************
+
+          IF (.NOT. ok_veget) THEN
+!
+! Read albedo
+             ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <ALB> est absent',1)
+
+             ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <ALB>',1)
+!
+! Read rugosity
+             ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <RUG> est absent',1)
+
+             ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
+             IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <RUG>',1)
+
+          END IF
+
+!****************************************************************************************
+! 5) Close file and distribuate variables to all processus
+!
+!****************************************************************************************
+          ierr = NF90_CLOSE(nid)
+          IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Pb when closing file', 1)
+       ENDIF ! is_mpi_root
+
+!$OMP END MASTER
+!$OMP BARRIER
+
+       IF ( type_ocean /= 'couple') THEN
+          CALL Scatter(sst_glo,sst)
+          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
+          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
+          IF (read_continents .OR. itime == 1) THEN
+             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
+             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
+          END IF
+       END IF
+
+       IF (.NOT. ok_veget) THEN
+          CALL Scatter(alb_glo, albedo)
+          CALL Scatter(rug_glo, rugos)
+       END IF
+
+    ENDIF ! time to read
+
+  END SUBROUTINE limit_read_tot
+
+
+END MODULE limit_read_mod
Index: LMDZ4/trunk/libf/phylmd/limit_slab.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/limit_slab.F90	(revision 996)
+++ LMDZ4/trunk/libf/phylmd/limit_slab.F90	(revision 996)
@@ -0,0 +1,121 @@
+! $Header$
+
+SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, lmt_foce, diff_sst)
+
+  USE dimphy
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_para
+  USE netcdf 
+
+  IMPLICIT NONE
+
+  INCLUDE "indicesol.h"
+  INCLUDE "temps.h"
+  INCLUDE "clesphys.h"
+  INCLUDE "dimensions.h"
+
+! In- and ouput arguments
+!****************************************************************************************
+  INTEGER, INTENT(IN) :: itime   ! numero du pas de temps courant
+  INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
+  REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
+  REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, lmt_foce, diff_sst
+
+! Locals variables with attribute SAVE
+!****************************************************************************************
+  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, foce_save
+
+! Locals variables
+!****************************************************************************************
+  INTEGER                  :: lmt_pas   
+  INTEGER                  :: nvarid, nid, ierr, i
+  INTEGER, DIMENSION(2)    :: start, epais 
+  REAL, DIMENSION(klon_glo):: bils_glo, foce_glo, sst_l_glo, sst_lp1_glo, diff_sst_glo
+  CHARACTER (len = 20)     :: modname = 'limit_slab'
+
+! End declaration
+!****************************************************************************************
+
+  ! calculate number of time steps for one day
+  lmt_pas = NINT(86400./dtime)
+  
+  IF (MOD(itime-1, lmt_pas) == 0) THEN   ! time to read
+     !$OMP MASTER  ! Only master thread
+     IF (is_mpi_root) THEN ! Only master processus
+        print*,'in limit_slab time to read, itime=',itime
+        
+        ierr = NF90_OPEN ('limit_slab.nc', NF90_NOWRITE, nid)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,&
+             'Pb in opening file limit_slab.nc',1)
+        
+        ! La tranche de donnees a lire:
+        start(1) = 1
+        start(2) = jour
+        epais(1) = klon_glo
+        epais(2) = 1
+
+!****************************************************************************************
+! 2) Read bils and ocean fraction
+!
+!****************************************************************************************
+!
+! Read bils_glo
+        ierr = NF90_INQ_VARID(nid, 'BILS_OCE', nvarid)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'The variable <BILS_OCE> is abstent',1)
+
+        ierr = NF90_GET_VAR(nid,nvarid,bils_glo,start,epais)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Reading of <BILS_OCE> failed',1)
+!
+! Read foce_glo
+        ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'The variable <FOCE> is abstent',1)
+
+        ierr = NF90_GET_VAR(nid,nvarid,foce_glo,start,epais)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Reading of <FOCE> failed',1)
+!
+! Read sst_glo for this day
+        ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'The variable <SST> is abstent',1)
+
+        ierr = NF90_GET_VAR(nid,nvarid,sst_l_glo,start,epais)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Reading of <SST> failed',1)
+
+! Read sst_glo for one day ahead
+        start(2) = jour + 1
+        IF (start(2) > 360) start(2)=1
+        ierr = NF90_GET_VAR(nid,nvarid,sst_lp1_glo,start,epais)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Reading of <SST> day+1 failed',1)
+
+! Calculate difference in temperature between this day and one ahead
+        DO i=1, klon_glo-1
+           diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
+        END DO
+        diff_sst_glo(klon_glo) = sst_lp1_glo(klon_glo) - sst_l_glo(1)
+
+!****************************************************************************************
+! 5) Close file and distribuate variables to all processus
+!
+!****************************************************************************************
+        ierr = NF90_CLOSE(nid)
+        IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Pb when closing file', 1)
+     ENDIF ! is_mpi_root
+
+!$OMP END MASTER
+       
+     IF (.NOT. ALLOCATED(bils_save)) THEN
+        ALLOCATE(bils_save(klon), foce_save(klon), stat=ierr)
+        IF (ierr /= 0) CALL abort_gcm('limit_slab', 'pb in allocation',1)
+     END IF
+
+     CALL Scatter(bils_glo, bils_save)
+     CALL Scatter(foce_glo, foce_save)
+     CALL Scatter(diff_sst_glo, diff_sst)
+     
+  ELSE ! not time to read
+     diff_sst(:) = 0.
+  ENDIF ! time to read
+
+  lmt_bils(:) = bils_save(:)
+  lmt_foce(:) = foce_save(:)
+  
+END SUBROUTINE limit_slab
Index: LMDZ4/trunk/libf/phylmd/oasis.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/oasis.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/oasis.F90	(revision 996)
@@ -44,5 +44,4 @@
   !$OMP THREADPRIVATE(out_var_id)
 
-  CHARACTER(LEN=*),PARAMETER :: OPA_version='OPA9'
 
 #ifdef CPP_COUPLE
@@ -58,4 +57,5 @@
 !     LF 09/2003
 !
+    USE surface_data, ONLY : version_ocean
     INCLUDE "dimensions.h"
 
@@ -130,14 +130,16 @@
 !     Define symbolic name for fields exchanged from atmos to coupler,
 !         must be the same as (1) of the field  definition in namcouple:
-    IF (OPA_version=='OPA9') THEN
-      cl_writ(1)='COTAUXXU'
-      cl_writ(2)='COTAUYYU'
-      cl_writ(3)='COTAUZZU'
-      cl_writ(4)='COTAUXXV'
-      cl_writ(5)='COTAUYYV'
-      cl_writ(6)='COTAUZZV'
-      cl_writ(7)='COWINDSP'
-      cl_writ(8)='COPEFWAT'
-      cl_writ(9)='COPEFICE'
+
+    cl_writ(1)='COTAUXXU'
+    cl_writ(2)='COTAUYYU'
+    cl_writ(3)='COTAUZZU'
+    cl_writ(4)='COTAUXXV'
+    cl_writ(5)='COTAUYYV'
+    cl_writ(6)='COTAUZZV'
+    cl_writ(7)='COWINDSP'
+    
+    IF (version_ocean=='nemo') THEN
+      cl_writ(8) ='COPEFWAT'
+      cl_writ(9) ='COPEFICE'
       cl_writ(10)='COTOSPSU'
       cl_writ(11)='COICEVAP'
@@ -150,14 +152,7 @@
       cl_writ(18)='CRWOCERD'
       cl_writ(19)='CRWOCECD'
-    ELSE IF (OPA_version=='OPA8') THEN
-      cl_writ(1)='COTAUXXU'
-      cl_writ(2)='COTAUYYU'
-      cl_writ(3)='COTAUZZU'
-      cl_writ(4)='COTAUXXV'
-      cl_writ(5)='COTAUYYV'
-      cl_writ(6)='COTAUZZV'
-      cl_writ(7)='COWINDSP'
-      cl_writ(8)='COSHFICE'
-      cl_writ(9)='COSHFOCE'
+    ELSE IF (version_ocean=='opa8') THEN
+      cl_writ(8) ='COSHFICE'
+      cl_writ(9) ='COSHFOCE'
       cl_writ(10)='CONSFICE'
       cl_writ(11)='CONSFOCE'
@@ -170,6 +165,4 @@
       cl_writ(18)='CORIVFLU'
       cl_writ(19)='COCALVIN'
-    ELSE
-      STOP 'Bad OPA version for coupled model'
     ENDIF
 
@@ -178,18 +171,16 @@
 !         must be the same as (2) of the field  definition in namcouple:
 !
-    IF (OPA_version=='OPA9') THEN
-      cl_read(1)='SISUTESW'
-      cl_read(2)='SIICECOV'
-      cl_read(4)='SIICEALW'
-      cl_read(3)='SIICTEMW'
-    ELSE IF (OPA_version=='OPA8') THEN
-      cl_read(1)='SISUTESW'
-      cl_read(2)='SIICECOV'
-      cl_read(3)='SIICEALW'
-      cl_read(4)='SIICTEMW'
-    ELSE
-      STOP 'Bad OPA version for coupled model'
-    ENDIF
-    
+    IF (version_ocean=='nemo') THEN
+       cl_read(1)='SISUTESW'
+       cl_read(2)='SIICECOV'
+       cl_read(4)='SIICEALW'
+       cl_read(3)='SIICTEMW'
+    ELSE IF (version_ocean=='opa8') THEN
+       cl_read(1)='SISUTESW'
+       cl_read(2)='SIICECOV'
+       cl_read(3)='SIICEALW'
+       cl_read(4)='SIICTEMW'
+    END IF
+
     il_var_nodims(1) = 2
     il_var_nodims(2) = 1
Index: LMDZ4/trunk/libf/phylmd/ocean_cpl_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/ocean_cpl_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/ocean_cpl_mod.F90	(revision 996)
@@ -16,10 +16,5 @@
   PRIVATE
 
-  PUBLIC :: ocean_cpl_init, ocean_cpl_get_vars, ocean_cpl_noice, ocean_cpl_ice
-
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE       :: tmp_flux_o
-  !$OMP THREADPRIVATE(tmp_flux_o)
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE       :: tmp_flux_g
-  !$OMP THREADPRIVATE(tmp_flux_g)
+  PUBLIC :: ocean_cpl_init, ocean_cpl_noice, ocean_cpl_ice
 
 !****************************************************************************************
@@ -43,17 +38,4 @@
     CHARACTER (len = 80) :: abort_message
     CHARACTER (len = 20) :: modname = 'ocean_cpl_init'
-
-
-    ALLOCATE(tmp_flux_o(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_o'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-
-    ALLOCATE(tmp_flux_g(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_g'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
 
 ! Initialize module cpl_init
@@ -71,8 +53,8 @@
        p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
        petAcoef, peqAcoef, petBcoef, peqBcoef, &
-       ps, u1_lay, v1_lay, pctsrf_in, &
+       ps, u1_lay, v1_lay, &
        radsol, snow, agesno, &
        qsurf, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+       tsurf_new, dflux_s, dflux_l)
 !
 ! This subroutine treats the "open ocean", all grid points that are not entierly covered
@@ -101,5 +83,4 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: ps
     REAL, DIMENSION(klon), INTENT(IN)        :: u1_lay, v1_lay
-    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf_in
 
 ! In/Output arguments
@@ -115,5 +96,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
 
 ! Local variables
@@ -122,5 +102,4 @@
     INTEGER, DIMENSION(1) :: iloc
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd
-    REAL, DIMENSION(klon) :: zx_sl
     REAL, DIMENSION(klon) :: fder_new
     REAL, DIMENSION(klon) :: tsurf_cpl
@@ -134,10 +113,8 @@
 
 !****************************************************************************************
-! Receive sea-surface temperature(tsurf_cpl) and new fraction of ocean surface(pctsrf_oce) 
-! from coupler
-!
-!****************************************************************************************
-    CALL cpl_receive_ocean_fields(itime, dtime, knon, knindex, pctsrf_in, &
-         tsurf_cpl, pctsrf_oce)
+! Receive sea-surface temperature(tsurf_cpl) from coupler
+!
+!****************************************************************************************
+    CALL cpl_receive_ocean_fields(knon, knindex, tsurf_cpl)
 
 !****************************************************************************************
@@ -176,21 +153,4 @@
 
 !****************************************************************************************
-! Flux ocean-atmosphere useful for "slab" ocean but here calculated only for printing
-! usage later in physiq  
-! 
-!****************************************************************************************
-    tmp_flux_o(:) = 0.0
-    DO i=1, knon
-       zx_sl(i) = RLVTT
-       IF (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
-       !IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
-       !       flux_o(i) = fluxsens(i) + fluxlat(i)
-       IF (pctsrf_oce(knindex(i)) .GT. epsfra) THEN
-          tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
-       ENDIF
-    ENDDO
-
-
-!****************************************************************************************
 ! Send and cumulate fields to the coupler
 !
@@ -213,8 +173,8 @@
        p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
        petAcoef, peqAcoef, petBcoef, peqBcoef, &
-       ps, u1_lay, v1_lay, pctsrf_in, &
+       ps, u1_lay, v1_lay, pctsrf, &
        radsol, snow, qsurf, &
        alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
+       tsurf_new, dflux_s, dflux_l)
 !
 ! This subroutine treats the ocean where there is ice. The subroutine first receives 
@@ -244,5 +204,5 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: ps
     REAL, DIMENSION(klon), INTENT(IN)        :: u1_lay, v1_lay
-    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf_in
+    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
 
 ! In/output arguments
@@ -258,5 +218,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_sic
 
 ! Local variables
@@ -277,11 +236,10 @@
 
 !****************************************************************************************
-! Receive ocean temperature(tsurf_cpl), albedo(alb_cpl) and new fraction of 
-! seaice(pctsrf_sic) from coupler
+! Receive ocean temperature(tsurf_cpl) and albedo(alb_new) from coupler
 !
 !****************************************************************************************
 
     CALL cpl_receive_seaice_fields(knon, knindex, &
-         tsurf_cpl, alb_cpl, pctsrf_sic)
+         tsurf_cpl, alb_cpl)
 
     alb1_new(1:knon) = alb_cpl(1:knon)
@@ -308,20 +266,5 @@
     CALL calcul_wind_flux(knon, dtime, taux, tauy)
     
-!****************************************************************************************
-! Flux ocean-atmosphere useful for "slab" ocean but here calculated only for printing
-! usage later in physiq  
-! 
-!  IM: faire dependre le coefficient de conduction de la glace de mer
-!      de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
-!      actuel correspond a 3m de glace de mer, cf. L.Li
-!
-!****************************************************************************************
-    tmp_flux_g(:) = 0.0
-    DO i = 1, knon
-       IF (cal(i) .GT. 1.0e-15 .AND. pctsrf_sic(knindex(i)) .GT. epsfra) &
-            tmp_flux_g(knindex(i)) = (tsurf_new(i)-t_grnd) * &
-            dif_grnd(i) * RCPD/cal(i) 
-    ENDDO
-    
+   
 !****************************************************************************************
 ! Calculate fder : flux derivative (sensible and latente)
@@ -344,5 +287,5 @@
 
     CALL cpl_send_seaice_fields(itime, dtime, knon, knindex, &
-       pctsrf_in, lafin, rlon, rlat, &
+       pctsrf, lafin, rlon, rlat, &
        swnet, lwnet, fluxlat, fluxsens, &
        precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, taux, tauy)
@@ -353,19 +296,3 @@
 !****************************************************************************************
 !
-  SUBROUTINE ocean_cpl_get_vars(flux_o, flux_g)
-
-! This subroutine returns variables private in this module to an external 
-! routine (physiq).
-
-    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o
-    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g
-
-! Set the output variables
-    flux_o(:) = tmp_flux_o(:)
-    flux_g(:) = tmp_flux_g(:)
-
-  END SUBROUTINE ocean_cpl_get_vars
-!
-!****************************************************************************************
-!
 END MODULE ocean_cpl_mod
Index: LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 996)
@@ -14,50 +14,5 @@
   IMPLICIT NONE
 
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE   :: tmp_flux_o, tmp_flux_g
-  !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
-
 CONTAINS
-!
-!****************************************************************************************
-! 
-  SUBROUTINE ocean_forced_init
-! Allocate fields needed for this module
-!    
-    INTEGER              :: error
-    CHARACTER (len = 80) :: abort_message
-    CHARACTER (len = 20) :: modname = 'ocean_forced_init'
-!****************************************************************************************
-
-    ALLOCATE(tmp_flux_o(1:klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_o'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-
-    ALLOCATE(tmp_flux_g(1:klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_g'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF    
-
-
-  END SUBROUTINE ocean_forced_init
-!
-!****************************************************************************************
-!****************************************************************************************
-! 
-  SUBROUTINE ocean_forced_final
-! Allocate fields needed for this module
-!    
-    INTEGER              :: error
-    CHARACTER (len = 80) :: abort_message
-    CHARACTER (len = 20) :: modname = 'ocean_forced_init'
-!****************************************************************************************
-
-    DEALLOCATE(tmp_flux_o)
-    DEALLOCATE(tmp_flux_g)
-
-
-  END SUBROUTINE ocean_forced_final
 !
 !****************************************************************************************
@@ -71,11 +26,12 @@
        radsol, snow, agesno, & 
        qsurf, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+       tsurf_new, dflux_s, dflux_l)
 !
 ! This subroutine treats the "open ocean", all grid points that are not entierly covered
 ! by ice.
-! The routine reads data from climatologie file and does some calculations at the 
+! The routine receives data from climatologie file limit.nc and does some calculations at the 
 ! surface. 
 !
+    USE limit_read_mod
     INCLUDE "indicesol.h"
     INCLUDE "YOMCST.h"
@@ -108,5 +64,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
-    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
 
 ! Local variables
@@ -116,5 +71,4 @@
     REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
     LOGICAL                     :: check=.FALSE.
-    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim
 
 !****************************************************************************************
@@ -125,13 +79,9 @@
 !****************************************************************************************
 ! 1)    
-! Read from climatologie file SST and fraction of sub-surfaces
-!
-!****************************************************************************************
-! Get from file tsurf_lim and pctsrf_lim
-    CALL interfoce_lim(itime, dtime, jour, & 
-         knon, knindex, &
-         debut, &
-         tsurf_lim, pctsrf_lim)
-    
+! Read sea-surface temperature from file limit.nc
+!
+!****************************************************************************************
+    CALL limit_read_sst(knon,knindex,tsurf_lim)
+
 !****************************************************************************************
 ! 2)
@@ -154,30 +104,4 @@
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 
-!****************************************************************************************
-! 3)
-! Calculate tmp_flux_o
-!
-!****************************************************************************************
-!IM: flux ocean-atmosphere utile pour le "slab" ocean
-! The flux are written to hist file
-    tmp_flux_o(:) = 0.0
-    DO i=1, knon
-       zx_sl(i) = RLVTT
-       IF (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
-
-!IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
-!       flux_o(i) = fluxsens(i) + fluxlat(i)
-       IF (pctsrf_lim(knindex(i),is_oce) .GT. epsfra) THEN
-          tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
-       ENDIF
-    ENDDO
-
-!****************************************************************************************
-! 4)
-! Return the new values for the ocean fraction
-!
-!****************************************************************************************
-
-    pctsrf_oce(:) = pctsrf_lim(:,is_oce)
 
   END SUBROUTINE ocean_forced_noice
@@ -192,10 +116,12 @@
        radsol, snow, qsol, agesno, tsoil, &
        qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
+       tsurf_new, dflux_s, dflux_l)
 !
 ! This subroutine treats the ocean where there is ice. 
 ! The routine reads data from climatologie file and does flux calculations at the 
 ! surface.
-!    
+!
+    USE limit_read_mod
+
     INCLUDE "indicesol.h"
     INCLUDE "dimsoil.h"
@@ -219,5 +145,4 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
 
-
 ! In/Output arguments
 !****************************************************************************************
@@ -226,5 +151,4 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
-
 
 ! Output arguments
@@ -236,6 +160,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic
-
 
 ! Local variables
@@ -248,5 +170,4 @@
     REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
     REAL, DIMENSION(klon)       :: soilcap, soilflux
-    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim
 
 !****************************************************************************************
@@ -256,35 +177,17 @@
 
 !****************************************************************************************
-! 1)    
-! Read from climatologie file SST and fraction of sub-surfaces
-!
-!****************************************************************************************
-    CALL interfoce_lim(itime, dtime, jour, & 
-         knon, knindex, &
-         debut, &
-         tsurf_tmp, pctsrf_lim)
-  
-    DO i = 1, knon
-       tsurf_tmp(i) = tsurf_in(i)
-       IF (pctsrf_lim(knindex(i),is_sic) < EPSFRA) THEN
-          snow(i) = 0.0
-          tsurf_tmp(i) = RTT - 1.8
-          IF (soil_model) tsoil(i,:) = RTT -1.8
-       ENDIF
-    ENDDO
-
-!****************************************************************************************
-! 2) 
+! 1) 
 ! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, 
 !                    dflux_s, dflux_l and qsurf
 !****************************************************************************************
+    tsurf_tmp(:) = tsurf_in(:)
 
 ! calculate the parameters cal, beta, capsol and dif_grnd
-
     CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
+
     
     IF (soil_model) THEN 
 ! update tsoil and calculate soilcap and soilflux
-       CALL soil(dtime, is_sic, knon,snow, tsurf_tmp, tsoil,soilcap, soilflux)
+       CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, tsoil,soilcap, soilflux)
        cal(1:knon) = RCPD / soilcap(1:knon)
        radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
@@ -305,5 +208,5 @@
 
 !****************************************************************************************
-! 3)
+! 2)
 ! Calculations due to snow and runoff
 !
@@ -327,49 +230,8 @@
     alb2_new(:) = alb1_new(:)
 
-!****************************************************************************************
-! 4)
-! Calculate tmp_flux_g
-!
-!****************************************************************************************
-!
-    tmp_flux_g(:) = 0.0
-    DO i = 1, knon
-!IM: faire dependre le coefficient de conduction de la glace de mer
-!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
-!    actuel correspond a 3m de glace de mer, cf. L.Li
-!
-       IF (cal(i) .GT. 1.0e-15 .AND. pctsrf_lim(knindex(i),is_sic) .GT. epsfra) &
-            tmp_flux_g(knindex(i)) = (tsurf_new(i)-t_grnd) * dif_grnd(i) *RCPD/cal(i)
-       
-    ENDDO
-
-
-!****************************************************************************************
-! 5)
-! Return the new values for the seaice fraction
-!
-!****************************************************************************************
-
-    pctsrf_sic(:) = pctsrf_lim(:,is_sic)
-
   END SUBROUTINE ocean_forced_ice
 !
 !****************************************************************************************
 !
-  SUBROUTINE ocean_forced_get_vars(flux_o, flux_g)
-! Get some variables from module oceanforced.
-! This subroutine returns variables to a external routine
-
-    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o
-    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g
-
-! Initialize the output variables
-    flux_o(:) = tmp_flux_o(:)
-    flux_g(:) = tmp_flux_g(:)
-
-  END SUBROUTINE ocean_forced_get_vars
-!
-!****************************************************************************************
-!
 END MODULE ocean_forced_mod
 
Index: LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 996)
@@ -7,5 +7,5 @@
 ! "ocean=slab".
 !
-  USE surface_data,     ONLY : tau_gl, calice, calsno
+  USE surface_data
   USE fonte_neige_mod,  ONLY : fonte_neige
   USE calcul_fluxs_mod, ONLY : calcul_fluxs
@@ -13,25 +13,6 @@
   
   IMPLICIT NONE
-
-  INTEGER, PRIVATE, SAVE                           :: lmt_pas, julien, idayvrai
-  !$OMP THREADPRIVATE(lmt_pas,julien,idayvrai)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_seaice
-  !$OMP THREADPRIVATE(tmp_seaice)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_tslab_loc
-  !$OMP THREADPRIVATE(tmp_tslab_loc)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slab_bils
-  !$OMP THREADPRIVATE(slab_bils)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE , SAVE  :: lmt_bils
-  !$OMP THREADPRIVATE(lmt_bils)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tmp_pctsrf_slab
-  !$OMP THREADPRIVATE(tmp_pctsrf_slab)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_tslab
-  !$OMP THREADPRIVATE(tmp_tslab)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_radsol
-  !$OMP THREADPRIVATE(tmp_radsol)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_flux_o, tmp_flux_g
-  !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
-  LOGICAL, PRIVATE, SAVE                           :: check = .FALSE.
-  !$OMP THREADPRIVATE(check)
+  PRIVATE
+  PUBLIC :: ocean_slab_frac, ocean_slab_noice
 
 CONTAINS
@@ -39,107 +20,45 @@
 !****************************************************************************************
 !
-  SUBROUTINE ocean_slab_init(dtime, tslab_rst, seaice_rst, pctsrf_rst)
+  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
 
+    USE limit_read_mod
     INCLUDE "indicesol.h"
-    INCLUDE "iniprint.h"
+!    INCLUDE "clesphys.h"
 
-! Input variables
+! Arguments
 !****************************************************************************************
-    REAL, INTENT(IN)                         :: dtime
-! Variables read from restart file
-    REAL, DIMENSION(klon), INTENT(IN)        :: tslab_rst          
-    REAL, DIMENSION(klon), INTENT(IN)        :: seaice_rst
-    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
-
+    INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
+    INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
+    REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf  ! sub-surface fraction
+    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
 
 ! Local variables
 !****************************************************************************************
-    INTEGER                :: error
     CHARACTER (len = 80)   :: abort_message
-    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
+    CHARACTER (len = 20)   :: modname = 'ocean_slab_frac'
 
 
-    WRITE(lunout,*) '************************'
-    WRITE(lunout,*) 'SLAB OCEAN est actif, prenez precautions !'
-    WRITE(lunout,*) '************************'    
+    IF (version_ocean=='sicOBS') THEN   
+       CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
+    ELSE
+       abort_message='Ocean slab model without forced sea-ice fractions has to be rewritten!!!'
+       CALL abort_gcm(modname,abort_message,1)
+! Here should sea-ice/ocean fraction either be calculated or returned if saved as a module varaiable 
+! (in the case the new fractions are calculated in ocean_slab_ice or ocean_slab_noice subroutines).  
+    END IF
 
-! Allocate variables initialize from restart fields
-    ALLOCATE(tmp_tslab(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_tslab'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-    tmp_tslab(:) = tslab_rst(:)
-
-    ALLOCATE(tmp_tslab_loc(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_tslab_loc'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-    tmp_tslab_loc(:) = tslab_rst(:)
-
-    ALLOCATE(tmp_seaice(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_seaice'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-    tmp_seaice(:) = seaice_rst(:)
-
-    ALLOCATE(tmp_pctsrf_slab(klon,nbsrf), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_pctsrf_slab'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-    tmp_pctsrf_slab(:,:) = pctsrf_rst(:,:)
-    
-! Allocate some other variables internal in module mod_oceanslab
-    ALLOCATE(tmp_radsol(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_radsol'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-
-    ALLOCATE(tmp_flux_o(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_o'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-    
-    ALLOCATE(tmp_flux_g(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_g'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-
-! a mettre un slab_bils aussi en force !!!
-    ALLOCATE(slab_bils(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation slab_bils'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-    slab_bils(:) = 0.0   
-
-    ALLOCATE(lmt_bils(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation lmt_bils'
-       CALL abort_gcm(modname,abort_message,1)
-    ENDIF
-
-
-! pour une lecture une fois par jour    
-    lmt_pas = NINT(86400./dtime * 1.0) 
-
-  END SUBROUTINE ocean_slab_init
+  END SUBROUTINE ocean_slab_frac
 !
 !****************************************************************************************
 !
   SUBROUTINE ocean_slab_noice( & 
-       dtime, knon, knindex, &
+       itime, dtime, jour, knon, knindex, &
        p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
        petAcoef, peqAcoef, petBcoef, peqBcoef, &
-       ps, u1_lay, v1_lay, &
+       ps, u1_lay, v1_lay, tsurf_in, &
        radsol, snow, agesno, &
        qsurf, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+       tsurf_new, dflux_s, dflux_l, lmt_bils)
 
     INCLUDE "indicesol.h"
@@ -148,4 +67,6 @@
 ! Input arguments
 !****************************************************************************************
+    INTEGER, INTENT(IN)                  :: itime
+    INTEGER, INTENT(IN)                  :: jour
     INTEGER, INTENT(IN)                  :: knon
     INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
@@ -159,4 +80,5 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: ps
     REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
+    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
 
 ! In/Output arguments
@@ -172,473 +94,54 @@
     REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)   :: pctsrf_oce
+    REAL, DIMENSION(klon), INTENT(OUT)   :: lmt_bils
 
 ! Local variables
 !****************************************************************************************
-    INTEGER                :: i
-    REAL, DIMENSION(klon)  :: cal, beta, dif_grnd
-    REAL, DIMENSION(klon)  :: tsurf_temp
+    INTEGER               :: i
+    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
+    REAL, DIMENSION(klon) :: lmt_bils_oce, lmt_foce, diff_sst
+    REAL                  :: calc_bils_oce, deltat
+    REAL, PARAMETER       :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
 
 !****************************************************************************************
-    IF (check) WRITE(*,*)' Entering ocean_slab_noice'    
-
-    tsurf_new(1:knon) = tmp_tslab(knindex(1:knon))
-    pctsrf_oce(:)   = tmp_pctsrf_slab(:,is_oce)
-    
-    tsurf_temp(:) = tsurf_new(:)
-    cal = 0.
-    beta = 1.
-    dif_grnd = 0.
-    agesno(:) = 0.
+! 1) Flux calculation
+!
+!****************************************************************************************
+    cal(:)      = 0.
+    beta(:)     = 1.
+    dif_grnd(:) = 0.
+    agesno(:)   = 0.
     
     CALL calcul_fluxs(knon, is_oce, dtime, &
-         tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
+         tsurf_in, p1lay, cal, beta, tq_cdrag, ps, &
          precip_rain, precip_snow, snow, qsurf,  &
          radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
          petAcoef, peqAcoef, petBcoef, peqBcoef, &
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
-    
-    tmp_flux_o(:) = 0.0
-    tmp_radsol(:) = 0.0
 
-    DO i=1, knon
-       tmp_radsol(knindex(i))=radsol(i)
-       
-       IF (pctsrf_oce(knindex(i)) .GT. epsfra) &
-            tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
-    ENDDO
-    
+!****************************************************************************************
+! 2) Get global variables lmt_bils and lmt_foce from file limit_slab.nc
+!
+!****************************************************************************************
+    CALL limit_slab(itime, dtime, jour, lmt_bils, lmt_foce, diff_sst)  ! global pour un processus
+
+    lmt_bils_oce(:) = 0.
+    WHERE (lmt_foce > 0.) 
+       lmt_bils_oce = lmt_bils / lmt_foce ! global 
+    END WHERE
+
+!****************************************************************************************
+! 3) Recalculate new temperature
+!
+!****************************************************************************************
+    DO i = 1, knon
+       calc_bils_oce = radsol(i) + fluxsens(i) + fluxlat(i)
+       deltat        = (calc_bils_oce - lmt_bils_oce(knindex(i)))*dtime/cyang +diff_sst(knindex(i))
+       tsurf_new(i)  = tsurf_in(i) + deltat
+    END DO
+
   END SUBROUTINE ocean_slab_noice
 !
 !****************************************************************************************
 !
-  SUBROUTINE ocean_slab_ice(   &
-       itime, dtime, jour, knon, knindex, &
-       debut, &
-       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
-       petAcoef, peqAcoef, petBcoef, peqBcoef, &
-       ps, u1_lay, v1_lay, &
-       radsol, snow, qsurf, qsol, agesno, tsoil, &
-       alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
-
-    INCLUDE "indicesol.h"
-    INCLUDE "dimsoil.h"
-    INCLUDE "YOMCST.h"
-    INCLUDE "iniprint.h"
-    INCLUDE "clesphys.h"
-
-! Input arguments  
-!****************************************************************************************
-    INTEGER, INTENT(IN)                  :: itime, jour, knon
-    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
-    REAL, INTENT(IN)                     :: dtime
-    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf
-    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
-    REAL, DIMENSION(klon), INTENT(IN)    :: tq_cdrag
-    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
-    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
-    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
-    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
-    REAL, DIMENSION(klon), INTENT(IN)    :: ps
-    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
-    LOGICAL, INTENT(IN)                  :: debut
-
-!In/Output arguments
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
-    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
-
-! Output arguments
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
-    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
-    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
-    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
-    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic
-
-! Local variables
-!****************************************************************************************
-    INTEGER                              :: i
-    REAL, DIMENSION(klon)                :: cal, beta, dif_grnd, capsol
-    REAL, DIMENSION(klon)                :: alb_neig, tsurf_temp
-    REAL, DIMENSION(klon)                :: soilcap, soilflux
-    REAL, DIMENSION(klon)                :: zfra
-    REAL, PARAMETER                      :: t_grnd=271.35
-    REAL                                 :: amn, amx 
-    REAL, DIMENSION(klon)                :: tslab
-    REAL, DIMENSION(klon)                :: seaice ! glace de mer (kg/m2)
-    REAL, DIMENSION(klon,nbsrf)          :: pctsrf_new
-
-!****************************************************************************************
-
-    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 
-
-! Initialization of output variables
-    alb1_new(:) = 0.0
-
-!****************************************************************************************
-! 
-!
-!****************************************************************************************
-    IF ( ok_slab_sicOBS) THEN   
-       ! glace de mer observee, lecture conditions limites
-       CALL interfoce_lim(itime, dtime, jour, & 
-            knon, knindex, &
-            debut, &
-            tsurf_new, pctsrf_new)
-
-       tmp_pctsrf_slab(:,:) = pctsrf_new(:,:)
-       WRITE(lunout,*) 'jour lecture pctsrf_new sic =',jour
-
-    ELSE
-       pctsrf_new=tmp_pctsrf_slab
-    ENDIF
-
-    DO i = 1, knon
-       tsurf_new(i) = tsurf(i)
-       IF (pctsrf_new(knindex(i),is_sic) < EPSFRA) THEN
-          snow(i) = 0.0
-          tsurf_new(i) = RTT - 1.8
-          IF (soil_model) tsoil(i,:) = RTT -1.8
-       ENDIF
-    ENDDO
-    
-    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
-    
-    IF (soil_model) THEN 
-       CALL soil(dtime, is_sic, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
-       cal(1:knon) = RCPD / soilcap(1:knon)
-       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
-    ELSE 
-       dif_grnd = 1.0 / tau_gl
-       cal = RCPD * calice
-       WHERE (snow > 0.0) cal = RCPD * calsno 
-    ENDIF
-    tsurf_temp = tsurf_new
-    beta = 1.0
-
-!****************************************************************************************
-!
-!
-!****************************************************************************************
-    CALL calcul_fluxs(knon, is_sic, dtime, &
-         tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
-         precip_rain, precip_snow, snow, qsurf,  &
-         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
-         petAcoef, peqAcoef, petBcoef, peqBcoef, &
-         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
-    
-    CALL fonte_neige( knon, is_sic, knindex, dtime, &
-         tsurf_temp, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap)
-
-!****************************************************************************************
-!     calcul albedo
-!
-!****************************************************************************************
-    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
-    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
-    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
-    alb1_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & 
-         0.6 * (1.0-zfra(1:knon))
-    
-    alb2_new(:) = alb1_new(:)
-
-!
-!IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
-    tmp_flux_g(:) = 0.0
-    DO i = 1, knon
-!
-!IM: faire dependre le coefficient de conduction de la glace de mer
-!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
-!    actuel correspond a 3m de glace de mer, cf. L.Li
-!
-       IF ((cal(i).GT.1.0e-15) .AND. (pctsrf_new(knindex(i),is_sic) .GT. epsfra)) THEN
-          tmp_flux_g(knindex(i))=(tsurf_new(i)-t_grnd) &
-               * dif_grnd(i) *RCPD/cal(i)
-       ENDIF
-!
-!IM: Attention: ne pas initialiser le tmp_radsol puisque c'est deja fait sur is_oce;
-!IM:            tmp_radsol doit etre le flux solaire qui arrive sur l'ocean 
-!IM:            et non pas celui qui arrive sur la glace de mer
-    ENDDO
-    
-!
-! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
-!
-
-    IF (check) THEN
-       amn=MIN(tmp_tslab(1),1000.)
-       amx=MAX(tmp_tslab(1),-1000.)
-       DO i=2, klon
-          amn=MIN(tmp_tslab(i),amn)
-          amx=MAX(tmp_tslab(i),amx)
-       ENDDO
-       
-       WRITE(lunout,*) ' debut avant interfoce_slab min max tmp_tslab',amn,amx
-    ENDIF !(check) THEN
-
-    tslab = tmp_tslab   
-
-    CALL interfoce_slab(klon, debut, itime, dtime, jour, &
-         tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
-         tslab, seaice, pctsrf_new)
-    
-    tmp_pctsrf_slab=pctsrf_new
-
-    DO i=1, knon
-       tmp_tslab(knindex(i))=tslab(knindex(i))
-    ENDDO
-     
-        
-!****************************************************************************************
-! Return the fraction of sea-ice
-! NB! jg : Peut-etre un probleme, faut-il prend aussi tmp_pctsrf_slab(:,is_oce)???
-!****************************************************************************************
-    pctsrf_sic(:) =  tmp_pctsrf_slab(:,is_sic)
-
-
-  END SUBROUTINE ocean_slab_ice
-! 
-!****************************************************************************************
-!
-  SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &
-       radsol, fluxo, fluxg, pctsrf, &
-       tslab, seaice, pctsrf_slab)
-!
-! Cette routine calcule la temperature d'un slab ocean, la glace de mer 
-! et les pourcentages de la maille couverte par l'ocean libre et/ou 
-! la glace de mer pour un "slab" ocean de 50m
-!
-! Conception: Laurent Li
-! Re-ecriture + adaptation LMDZ4: I. Musat
-!
-! input:
-!   klon         nombre total de points de grille
-!   debut        logical: 1er appel a la physique
-!   itap         numero du pas de temps
-!   dtime        pas de temps de la physique (en s)
-!   ijour        jour dans l'annee en cours
-!   radsol       rayonnement net au sol (LW + SW)
-!   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques 
-!   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean
-!   pctsrf       tableau des pourcentages de surface de chaque maille
-! output: 
-!   tslab        temperature de l'ocean libre
-!   seaice       glace de mer (kg/m2)
-!   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
-
-    INCLUDE "indicesol.h"
-    INCLUDE "YOMCST.h"
-    INCLUDE "iniprint.h"
-    INCLUDE "clesphys.h"
-
-! Input arguments
-!****************************************************************************************
-    INTEGER, INTENT(IN)                       :: klon
-    LOGICAL, INTENT(IN)                       :: debut    ! not used
-    INTEGER, INTENT(IN)                       :: itap
-    REAL, INTENT(IN)                          :: dtime       ! not used
-    INTEGER, INTENT(IN)                       :: ijour
-    REAL, DIMENSION(klon), INTENT(IN)         :: radsol
-    REAL, DIMENSION(klon), INTENT(IN)         :: fluxo
-    REAL, DIMENSION(klon), INTENT(IN)         :: fluxg
-    REAL, DIMENSION(klon, nbsrf), INTENT(IN)  :: pctsrf
-
-! Output arguments
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(OUT)        :: tslab
-    REAL, DIMENSION(klon), INTENT(OUT)        :: seaice ! glace de mer (kg/m2)
-    REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: pctsrf_slab
-
-! Local variables
-!****************************************************************************************
-    REAL                    :: amn, amx
-    REAL, PARAMETER         :: unjour=86400.
-    REAL, PARAMETER         :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
-    REAL, PARAMETER         :: cbing=0.334e+05        ! J/kg
-    REAL, DIMENSION(klon)   :: siceh !hauteur de la glace de mer (m)
-    INTEGER                 :: i
-    REAL                    :: zz, za, zb
-!
-!****************************************************************************************
-!
-    julien = MOD(ijour,360)
-
-    IF (check ) THEN
-       amn=MIN(tmp_tslab_loc(1),1000.)
-       amx=MAX(tmp_tslab_loc(1),-1000.)
-       DO i=2, klon
-          amn=MIN(tmp_tslab_loc(i),amn)
-          amx=MAX(tmp_tslab_loc(i),amx)
-       ENDDO
-
-       WRITE(lunout,*) ' debut min max tslab',amn,amx 
-       WRITE(lunout,*) ' itap,lmt_pas unjour',itap,lmt_pas,unjour
-    ENDIF
-
-    pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)
-!
-! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
-! 
-    IF (MOD(itap,lmt_pas) .EQ. 1) THEN 
-       ! 1er pas de temps de la journee
-       idayvrai = ijour
-       CALL condsurf(julien,idayvrai, lmt_bils)
-    ENDIF
-
-    DO i = 1, klon
-       IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
-            (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
-!
-! fabriquer de la glace si congelation atteinte:
-!          
-          IF (tmp_tslab_loc(i).LT.(RTT-1.8)) THEN
-             zz =  (RTT-1.8)-tmp_tslab_loc(i)
-             tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
-             seaice(i) = tmp_seaice(i)
-             tmp_tslab_loc(i) = RTT-1.8
-          ENDIF
-!
-! faire fondre de la glace si temperature est superieure a 0:
-!
-          IF ((tmp_tslab_loc(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN
-             zz = cyang/cbing * (tmp_tslab_loc(i)-RTT)
-             zz = MIN(zz,tmp_seaice(i))
-             tmp_seaice(i) = tmp_seaice(i) - zz
-             seaice(i) = tmp_seaice(i)
-             tmp_tslab_loc(i) = tmp_tslab_loc(i) - zz*cbing/cyang
-          ENDIF
-!
-! limiter la glace de mer a 10 metres (10000 kg/m2)
-!
-          IF(tmp_seaice(i).GT.45.) THEN
-             tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)
-          ELSE
-             tmp_seaice(i) = 0. 
-          ENDIF
-          seaice(i) = tmp_seaice(i)
-          siceh(i)=tmp_seaice(i)/1000. !en metres
-!
-! determiner la nature du sol (glace de mer ou ocean libre):
-!
-! on fait dependre la fraction de seaice "pctsrf(i,is_sic)" 
-! de l'epaisseur de seaice :
-! pctsrf(i,is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm
-! et pctsrf(i,is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur
-!
-
-          IF(.NOT.ok_slab_sicOBS) THEN
-             pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &
-                  1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))
-             pctsrf_slab(i,is_oce)=1.0 - &
-                  (pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)+pctsrf_slab(i,is_sic))
-          ELSE
-             IF (i.EQ.1) WRITE(lunout,*) 'cas ok_slab_sicOBS TRUE : passe sur la modif.'
-          ENDIF !(.NOT.ok_slab_sicOBS) then
-       ENDIF !pctsrf
-    ENDDO
-!
-! Calculer le bilan du flux de chaleur au sol :
-!
-    DO i = 1, klon
-       za = radsol(i) + fluxo(i)
-       zb = fluxg(i)
-       IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
-            (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
-          slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &
-               +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)
-       ENDIF
-    ENDDO !klon
-!
-! calcul tslab 
-!
-    IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee
-!
-! calcul tslab 
-       amn=MIN(tmp_tslab_loc(1),1000.)
-       amx=MAX(tmp_tslab_loc(1),-1000.)
-       DO i = 1, klon
-          IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
-               (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
-             tmp_tslab_loc(i) = tmp_tslab_loc(i) + &
-                  (slab_bils(i)-lmt_bils(i)) &
-                  /cyang*unjour
-
-! on remet l'accumulation a 0
-             slab_bils(i) = 0.
-          ENDIF !pctsrf
-!
-          IF (check) THEN
-             IF(i.EQ.1) THEN  
-                WRITE(lunout,*) 'i tmp_tslab_loc MOD(itap,lmt_pas).EQ.0 sicINTER',i,itap, &
-                     tmp_tslab_loc(i), tmp_seaice(i)
-             ENDIF
-             
-             amn=MIN(tmp_tslab_loc(i),amn)
-             amx=MAX(tmp_tslab_loc(i),amx)
-          ENDIF
-       ENDDO !klon
-    ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN
-
-    IF ( check ) THEN
-       WRITE(lunout,*) 'fin min max tslab',amn,amx 
-    ENDIF
-
-    tslab  = tmp_tslab_loc
-    seaice = tmp_seaice
-
-  END SUBROUTINE interfoce_slab
-!
-!****************************************************************************************  
-!
-  SUBROUTINE ocean_slab_final(tslab_rst, seaice_rst)
-
-! This subroutine will send to phyredem the variables concerning the slab 
-! ocean that should be written to restart file.
-
-!****************************************************************************************
-
-    REAL, DIMENSION(klon), INTENT(OUT) :: tslab_rst
-    REAL, DIMENSION(klon), INTENT(OUT) :: seaice_rst
-
-!****************************************************************************************
-! Set the output variables
-    tslab_rst(:)  = tmp_tslab(:)
-!    tslab_rst(:)  = tmp_tslab_loc(:)
-    seaice_rst(:) = tmp_seaice(:)
-
-! Deallocation of all variables in module
-    DEALLOCATE(tmp_tslab, tmp_tslab_loc, tmp_pctsrf_slab)
-    DEALLOCATE(tmp_seaice, tmp_radsol, tmp_flux_o, tmp_flux_g)
-    DEALLOCATE(slab_bils, lmt_bils)
-
-  END SUBROUTINE ocean_slab_final
-!
-!****************************************************************************************
-!
-  SUBROUTINE ocean_slab_get_vars(tslab_loc, seaice_loc, flux_o_loc, flux_g_loc)
-! "Get some variables from module ocean_slab_mod"
-! This subroutine prints variables to a external routine
-
-    REAL, DIMENSION(klon), INTENT(OUT) :: tslab_loc
-    REAL, DIMENSION(klon), INTENT(OUT) :: seaice_loc
-    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o_loc
-    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g_loc
-
-! Set the output variables
-    tslab_loc(:)  = tmp_tslab(:)
-    seaice_loc(:) = tmp_seaice(:)
-    flux_o_loc(:) = tmp_flux_o(:)
-    flux_g_loc(:) = tmp_flux_g(:)
-
-  END SUBROUTINE ocean_slab_get_vars
-!
-!****************************************************************************************
-!
 END MODULE ocean_slab_mod
Index: LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 996)
@@ -13,5 +13,5 @@
   USE mod_phys_lmdz_para,  ONLY : mpi_size
   USE ioipsl
-  USE surface_data,        ONLY : ocean, ok_veget
+  USE surface_data,        ONLY : type_ocean, ok_veget
   USE surf_land_mod,       ONLY : surf_land
   USE surf_landice_mod,    ONLY : surf_landice
@@ -151,7 +151,7 @@
 !****************************************************************************************
 
-    IF (ocean /= 'slab  ' .AND. ocean /= 'force ' .AND. ocean /= 'couple') THEN
+    IF (type_ocean /= 'slab  ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
       WRITE(lunout,*)' *** Warning ***'
-      WRITE(lunout,*)'Option couplage pour l''ocean = ', ocean
+      WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
       abort_message='option pour l''ocean non valable'
       CALL abort_gcm(modname,abort_message,1)
@@ -187,5 +187,5 @@
        zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
        d_t,       d_q,       d_u,      d_v,           & 
-       zcoefh,    pctsrf_new,                         &
+       zcoefh,    slab_wfbils,                        &
        qsol_d,    zq2m,      s_pblh,   s_plcl,        &
        s_capCL,   s_oliqCL,  s_cteiCL, s_pblT,        &
@@ -299,4 +299,5 @@
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
+    REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke 
 
 ! Output variables
@@ -321,7 +322,7 @@
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zcoefh     ! coef for turbulent diffusion of T and Q, mean for each grid point
-    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pctsrf_new ! new sub-surface fraction
 
 ! Output only for diagnostics
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: slab_wfbils! heat balance at surface only for slab at ocean points
     REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol_d     ! water height in the soil (mm)
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m       ! water vapour at 2m, mean for each grid point
@@ -367,7 +368,4 @@
     REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m        ! water vapour at 2 meter height
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
-
-! Input/output
-    REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke
 
 
@@ -431,5 +429,4 @@
     REAL, DIMENSION(klon)              :: ypsref
     REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb1_new, yalb2_new
-    REAL, DIMENSION(klon)              :: pctsrf_nsrf
     REAL, DIMENSION(klon)              :: ztsol
     REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
@@ -446,5 +443,4 @@
     REAL, DIMENSION(klon,klev+1)       :: ytke
     REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
-    REAL, DIMENSION(klon,nbsrf)        :: pctsrf_pot
     CHARACTER(len=80)                  :: abort_message
     CHARACTER(len=20)                  :: modname = 'pbl_surface'
@@ -470,16 +466,10 @@
     REAL, DIMENSION(klon,nbsrf)        :: trmb2        ! inhibition
     REAL, DIMENSION(klon,nbsrf)        :: trmb3        ! point Omega
-    REAL, DIMENSION(klon,nbsrf)        :: zx_rh2m, zx_qsat2m
-    REAL, DIMENSION(klon,nbsrf)        :: zx_qs1, zx_t1
-    REAL, DIMENSION(klon,nbsrf)        :: zdelta1, zcor1
+    REAL, DIMENSION(klon,nbsrf)        :: zx_t1
     REAL, DIMENSION(klon, nbsrf)       :: alb          ! mean albedo for whole SW interval
     REAL, DIMENSION(klon)              :: ylwdown      ! jg : temporary (ysollwdown)
 
-
-!jg+ temporary test
-    REAL, DIMENSION(klon,klev)         :: y_flux_u_old, y_flux_v_old
-    REAL, DIMENSION(klon,klev)         :: y_d_u_old, y_d_v_old
-!jg-
-    
+    REAL                               :: zx_qs1, zcor1, zdelta1 
+
 !****************************************************************************************
 ! Declarations specifiques pour le 1D. A reprendre 
@@ -558,5 +548,5 @@
     yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
     ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
-    yq = 0.0      ; pctsrf_new = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0 
+    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0 
     yrugoro = 0.0 ; yu10mx = 0.0     ; yu10my = 0.0    ; ywindsp = 0.0   
     d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
@@ -661,14 +651,8 @@
 ! 4) Loop over different surfaces
 !
-! All points with a possibility of the current surface are used. This is done
-! to allow the sea-ice to appear or disappear. It is considered here that the 
-! entier domaine of the ocean possibly can contain sea-ice.
+! Only points containing a fraction of the sub surface will be threated.
 ! 
 !****************************************************************************************
-
-    pctsrf_pot = pctsrf
-    pctsrf_pot(:,is_oce) = 1. - zmasq(:)
-    pctsrf_pot(:,is_sic) = 1. - zmasq(:)
-      
+    
     loop_nbsrf: DO nsrf = 1, nbsrf
 
@@ -677,5 +661,5 @@
        knon  = 0
        DO i = 1, klon
-          IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
+          IF (pctsrf(i,nsrf) > 0.) THEN
              knon = knon + 1
              ni(knon) = i
@@ -730,6 +714,6 @@
              ypaprs(j,k) = paprs(i,k)
              ypplay(j,k) = pplay(i,k)
-             ydelp(j,k) = delp(i,k)
-             ytke(j,k)=tke(i,k,nsrf)
+             ydelp(j,k)  = delp(i,k)
+             ytke(j,k)   = tke(i,k,nsrf)
              yu(j,k) = u(i,k)
              yv(j,k) = v(i,k)
@@ -762,5 +746,5 @@
        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf,  &
-            ycoefm, ycoefh,ytke)
+            ycoefm, ycoefh, ytke)
        
 !****************************************************************************************
@@ -817,5 +801,5 @@
                ysnow, yqsol, yagesno, ytsoil, &
                yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
-               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf, &
+               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
                ylwdown)
      
@@ -828,9 +812,9 @@
                ysnow, yqsurf, yqsol, yagesno, &
                ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
-               ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+               ytsurf_new, y_dflux_t, y_dflux_q)
           
        CASE(is_oce)
           CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
-               yrugos, ywindsp, rmu0, yfder, &
+               yrugos, ywindsp, rmu0, yfder, yts, &
                itap, dtime, jour, knon, ni, &
                debut, &
@@ -840,5 +824,5 @@
                ysnow, yqsurf, yagesno, &
                yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
-               ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+               ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils)
           
        CASE(is_sic)
@@ -852,5 +836,5 @@
                ysnow, yqsurf, yqsol, yagesno, ytsoil, &
                yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
-               ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+               ytsurf_new, y_dflux_t, y_dflux_q)
           
 
@@ -861,9 +845,4 @@
        END SELECT
 
-!****************************************************************************************
-! Save the fraction of this sub-surface 
-!
-!****************************************************************************************
-       pctsrf_new(:,nsrf) = pctsrf_nsrf(:)
 
 !****************************************************************************************
@@ -882,6 +861,6 @@
 !****************************************************************************************
 ! H and Q
-!      print *,'pbl_surface: ok_flux_surf=',ok_flux_surf
-!      print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
+!       print *,'pbl_surface: ok_flux_surf=',ok_flux_surf
+!       print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
        if (ok_flux_surf) then
           y_flux_t1(:) =  fsens
@@ -916,5 +895,5 @@
 !****************************************************************************************
 
-       tke(:,:,nsrf)=0.
+       tke(:,:,nsrf) = 0.
        DO k = 1, klev
           DO j = 1, knon
@@ -922,8 +901,8 @@
              ycoefh(j,k) = ycoefh(j,k) * ypct(j)
              ycoefm(j,k) = ycoefm(j,k) * ypct(j)
-             y_d_t(j,k) = y_d_t(j,k) * ypct(j)
-             y_d_q(j,k) = y_d_q(j,k) * ypct(j)
-             y_d_u(j,k) = y_d_u(j,k) * ypct(j)
-             y_d_v(j,k) = y_d_v(j,k) * ypct(j)
+             y_d_t(j,k)  = y_d_t(j,k) * ypct(j)
+             y_d_q(j,k)  = y_d_q(j,k) * ypct(j)
+             y_d_u(j,k)  = y_d_u(j,k) * ypct(j)
+             y_d_v(j,k)  = y_d_v(j,k) * ypct(j)
 
              flux_t(i,k,nsrf) = y_flux_t(j,k)
@@ -932,5 +911,5 @@
              flux_v(i,k,nsrf) = y_flux_v(j,k)
 
-             tke(i,k,nsrf)=ytke(j,k)
+             tke(i,k,nsrf)    = ytke(j,k)
 
           ENDDO
@@ -1021,9 +1000,11 @@
        trmb2(:,nsrf)  = 0.        ! inhibition 
        trmb3(:,nsrf)  = 0.        ! Point Omega
+       rh2m(:)        = 0.
+       qsat2m(:)      = 0.
 
 #undef T2m     
 #define T2m     
 #ifdef T2m
-! diagnostic t,q a 2m et u, v a 10m
+! Calculations of diagnostic t,q at 2m and u, v at 10m
 
        DO j=1, knon
@@ -1055,13 +1036,26 @@
           i = ni(j)
           t2m(i,nsrf)=yt2m(j)
+          q2m(i,nsrf)=yq2m(j)
           
-          q2m(i,nsrf)=yq2m(j)
-
-! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
+          ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
           u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
           v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
-
        END DO
 
+!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
+!IM Ajoute dependance type surface
+       IF (thermcep) THEN
+          DO j = 1, knon
+             i=ni(j)
+             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
+             zx_qs1  = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
+             zx_qs1  = MIN(0.5,zx_qs1)
+             zcor1   = 1./(1.-RETV*zx_qs1)
+             zx_qs1  = zx_qs1*zcor1
+             
+             rh2m(i)   = rh2m(i)   + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
+             qsat2m(i) = qsat2m(i) + zx_qs1  * pctsrf(i,nsrf)
+          END DO
+       END IF
 
        CALL HBTM(knon, ypaprs, ypplay, &
@@ -1086,27 +1080,7 @@
        
 #else 
-! not defined T2m
+! T2m not defined
 ! No calculation
-! Set output variables to zero
-       DO j = 1, knon
-          i = ni(j)
-          pblh(i,nsrf)   = 0.
-          plcl(i,nsrf)   = 0.
-          capCL(i,nsrf)  = 0.
-          oliqCL(i,nsrf) = 0.
-          cteiCL(i,nsrf) = 0.
-          pblT(i,nsrf)   = 0.
-          therm(i,nsrf)  = 0.
-          trmb1(i,nsrf)  = 0.
-          trmb2(i,nsrf)  = 0.
-          trmb3(i,nsrf)  = 0.
-       END DO
-       DO j = 1, knon
-          i = ni(j) 
-          t2m(i,nsrf)=0.
-          q2m(i,nsrf)=0.
-          u10m(i,nsrf)=0.
-          v10m(i,nsrf)=0.
-       END DO
+       PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
 #endif
 
@@ -1120,6 +1094,4 @@
 ! 16) Calculate the mean value over all sub-surfaces for som variables
 !
-!     NB!!! jg : Pour garder la convergence numerique j'utilise pctsrf_new comme c'etait 
-!     fait dans physiq.F mais ca devrait plutot etre pctsrf???!!!!! A verifier!
 !****************************************************************************************
     
@@ -1129,8 +1101,8 @@
        DO k = 1, klev
           DO i = 1, klon
-             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf_new(i,nsrf)
-             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf_new(i,nsrf)
-             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf_new(i,nsrf)
-             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf_new(i,nsrf)
+             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
+             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
+             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
+             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
           END DO
        END DO
@@ -1143,14 +1115,4 @@
     ENDDO
    
-
-    DO i = 1, klon
-       IF ( ABS( pctsrf_new(i, is_ter) + pctsrf_new(i, is_lic) + &
-            pctsrf_new(i, is_oce) + pctsrf_new(i, is_sic)  - 1.) .GT. EPSFRA) & 
-            THEN 
-          WRITE(*,*) 'physiq : pb sous surface au point ', i, &
-               pctsrf_new(i, 1 : nbsrf)
-       ENDIF
-    ENDDO
-
 !
 ! Incrementer la temperature du sol
@@ -1171,26 +1133,26 @@
           
           wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
-               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf_new(i,nsrf)
+               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
           wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
-               pctsrf_new(i,nsrf)
-
-          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf_new(i,nsrf)
-          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf_new(i,nsrf)
+               pctsrf(i,nsrf)
+
+          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf(i,nsrf)
+          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
           
-          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf_new(i,nsrf)
-          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf_new(i,nsrf)
-          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf_new(i,nsrf)
-          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf_new(i,nsrf)
-
-          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf_new(i,nsrf)
-          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf_new(i,nsrf)
-          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf_new(i,nsrf)
-          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf_new(i,nsrf)
-          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf_new(i,nsrf)
-          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf_new(i,nsrf)
-          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf_new(i,nsrf)
-          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf_new(i,nsrf)
-          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf_new(i,nsrf)
-          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf_new(i,nsrf)
+          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
+          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
+          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
+          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
+
+          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf(i,nsrf)
+          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf(i,nsrf)
+          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf(i,nsrf)
+          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
+          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
+          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf(i,nsrf)
+          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf(i,nsrf)
+          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf(i,nsrf)
+          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf(i,nsrf)
+          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf(i,nsrf)
        END DO
     END DO
@@ -1205,31 +1167,31 @@
        PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
     ENDIF
-!
-! If a sub-surface does not exsist for a grid point, the mean value for all 
-! sub-surfaces is distributed.
-!
-    DO nsrf = 1, nbsrf
-       DO i = 1, klon
-          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
-             ts(i,nsrf)     = zxtsol(i)
-             t2m(i,nsrf)    = zt2m(i)
-             q2m(i,nsrf)    = zq2m(i)
-             u10m(i,nsrf)   = zu10m(i)
-             v10m(i,nsrf)   = zv10m(i)
-
-! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
-             pblh(i,nsrf)   = s_pblh(i)
-             plcl(i,nsrf)   = s_plcl(i)
-             capCL(i,nsrf)  = s_capCL(i)
-             oliqCL(i,nsrf) = s_oliqCL(i) 
-             cteiCL(i,nsrf) = s_cteiCL(i)
-             pblT(i,nsrf)   = s_pblT(i)
-             therm(i,nsrf)  = s_therm(i)
-             trmb1(i,nsrf)  = s_trmb1(i)
-             trmb2(i,nsrf)  = s_trmb2(i)
-             trmb3(i,nsrf)  = s_trmb3(i)
-          ENDIF
-       ENDDO
-    ENDDO
+!!$!
+!!$! If a sub-surface does not exsist for a grid point, the mean value for all 
+!!$! sub-surfaces is distributed.
+!!$!
+!!$    DO nsrf = 1, nbsrf
+!!$       DO i = 1, klon
+!!$          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
+!!$             ts(i,nsrf)     = zxtsol(i)
+!!$             t2m(i,nsrf)    = zt2m(i)
+!!$             q2m(i,nsrf)    = zq2m(i)
+!!$             u10m(i,nsrf)   = zu10m(i)
+!!$             v10m(i,nsrf)   = zv10m(i)
+!!$
+!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
+!!$             pblh(i,nsrf)   = s_pblh(i)
+!!$             plcl(i,nsrf)   = s_plcl(i)
+!!$             capCL(i,nsrf)  = s_capCL(i)
+!!$             oliqCL(i,nsrf) = s_oliqCL(i) 
+!!$             cteiCL(i,nsrf) = s_cteiCL(i)
+!!$             pblT(i,nsrf)   = s_pblT(i)
+!!$             therm(i,nsrf)  = s_therm(i)
+!!$             trmb1(i,nsrf)  = s_trmb1(i)
+!!$             trmb2(i,nsrf)  = s_trmb2(i)
+!!$             trmb3(i,nsrf)  = s_trmb3(i)
+!!$          ENDIF
+!!$       ENDDO
+!!$    ENDDO
 
 
@@ -1242,32 +1204,9 @@
     DO nsrf = 1, nbsrf
        DO i = 1, klon
-          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf_new(i,nsrf)
-          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf_new(i,nsrf)
+          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
+          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
        END DO
     END DO
 
-!
-!IM Calculer l'humidite relative a 2m (rh2m) pour diagnostique
-!IM ajout dependance type surface
-    rh2m(:)   = 0.0
-    qsat2m(:) = 0.0
-
-    DO i = 1, klon
-       DO nsrf=1, nbsrf
-          zx_t1(i,nsrf) = t2m(i,nsrf)
-          IF (thermcep) THEN
-             zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf)))
-             zx_qs1(i,nsrf)  = r2es * &
-                  FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1)
-             zx_qs1(i,nsrf)  = MIN(0.5,zx_qs1(i,nsrf))
-             zcor1(i,nsrf)   = 1./(1.-retv*zx_qs1(i,nsrf))
-             zx_qs1(i,nsrf)  = zx_qs1(i,nsrf)*zcor1(i,nsrf)
-          END IF
-          zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf)
-          zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf)
-          rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf_new(i,nsrf)
-          qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf_new(i,nsrf)
-       END DO
-    END DO
 
 ! Some of the module declared variables are returned for printing in physiq.F
@@ -1322,5 +1261,111 @@
 !  
 !****************************************************************************************
+! 
+  SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, tke)
+
+    ! Give default values where new fraction has appread
+
+    INCLUDE "indicesol.h"
+    INCLUDE "dimsoil.h"
+    INCLUDE "clesphys.h"
+
+! Input variables
+!****************************************************************************************
+    INTEGER, INTENT(IN)                     :: itime
+    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
+
+! InOutput variables
+!****************************************************************************************
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: tsurf
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: alb1, alb2
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: u10m, v10m
+    REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
+
+! Local variables
+!****************************************************************************************
+    INTEGER           :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
+    CHARACTER(len=80) :: abort_message
+    CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
+    INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
+    LOGICAL           :: debug=.FALSE.
+!
+! All at once !! 
+!****************************************************************************************
+    
+    DO nsrf = 1, nbsrf
+       ! First decide complement sub-surfaces
+       SELECT CASE (nsrf)
+       CASE(is_oce)
+          nsrf_comp1=is_sic
+          nsrf_comp2=is_ter
+          nsrf_comp3=is_lic
+       CASE(is_sic)
+          nsrf_comp1=is_oce
+          nsrf_comp2=is_ter
+          nsrf_comp3=is_lic
+       CASE(is_ter)
+          nsrf_comp1=is_lic
+          nsrf_comp2=is_oce
+          nsrf_comp3=is_sic
+       CASE(is_lic)
+          nsrf_comp1=is_ter
+          nsrf_comp2=is_oce
+          nsrf_comp3=is_sic
+       END SELECT
+
+       ! Initialize all new fractions
+       DO i=1, klon
+          IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
+
+             IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
+                ! Use the complement sub-surface, keeping the continents unchanged
+                qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
+                evap(i,nsrf)  = evap(i,nsrf_comp1)
+                rugos(i,nsrf) = rugos(i,nsrf_comp1)
+                tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
+                alb1(i,nsrf)  = alb1(i,nsrf_comp1)
+                alb2(i,nsrf)  = alb2(i,nsrf_comp1)
+                u10m(i,nsrf)  = u10m(i,nsrf_comp1)
+                v10m(i,nsrf)  = v10m(i,nsrf_comp1)
+                tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
+                mfois(nsrf) = mfois(nsrf) + 1
+             ELSE
+                ! The continents have changed. The new fraction receives the mean sum of the existent fractions
+                qsurf(i,nsrf) = qsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
+                evap(i,nsrf)  = evap(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
+                rugos(i,nsrf) = rugos(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
+                tsurf(i,nsrf) = tsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
+                alb1(i,nsrf)  = alb1(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
+                alb2(i,nsrf)  = alb2(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
+                u10m(i,nsrf)  = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
+                v10m(i,nsrf)  = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
+                tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
+            
+                ! Security abort. This option has never been tested. To test, comment the following line.
+!                abort_message='The fraction of the continents have changed!'
+!                CALL abort_gcm(modname,abort_message,1)
+                nfois(nsrf) = nfois(nsrf) + 1
+             END IF
+             snow(i,nsrf)     = 0.
+             agesno(i,nsrf)   = 0.
+             ftsoil(i,:,nsrf) = tsurf(i,nsrf)
+          ELSE
+             pfois(nsrf) = pfois(nsrf)+ 1
+          END IF
+       END DO
+       
+    END DO
+
+    IF (debug) THEN
+       print*,'itime=,',itime, 'Pas de nouveau fraction',pfois,'fois'
+       print*,'itime=,',itime, 'The fraction of the continents have changed',nfois,'fois'
+       print*,'itime=,',itime, 'The fraction ocean-seaice has changed',mfois,'fois'
+    END IF
+
+  END SUBROUTINE pbl_surface_newfrac
+
 !  
+!****************************************************************************************
+!  
 
 END MODULE pbl_surface_mod
Index: LMDZ4/trunk/libf/phylmd/phyetat0.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/phyetat0.F	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/phyetat0.F	(revision 996)
@@ -5,5 +5,4 @@
 c
       SUBROUTINE phyetat0 (fichnom,
-     .           ocean_in, ok_veget_in,
      .           clesphy0,
      .           tabcntr0)
@@ -13,10 +12,8 @@
       USE mod_phys_lmdz_para
       USE iophy
-      USE ocean_slab_mod,   ONLY : ocean_slab_init
       USE ocean_cpl_mod,    ONLY : ocean_cpl_init
-      USE ocean_forced_mod, ONLY : ocean_forced_init
       USE fonte_neige_mod,  ONLY : fonte_neige_init
       USE pbl_surface_mod,  ONLY : pbl_surface_init
-      USE surface_data,     ONLY : ocean, ok_veget
+      USE surface_data,     ONLY : type_ocean
       USE phys_state_var_mod
 
@@ -68,5 +65,4 @@
       REAL wake_fip_glo(klon_glo)
       REAL tsoil_p(klon,nsoilmx,nbsrf)
-      REAL tslab_p(klon), seaice_p(klon)
       REAL qsurf_p(klon,nbsrf)
       REAL qsol_p(klon)
@@ -83,6 +79,4 @@
       REAL zmasq_glo(klon_glo)
       REAL tsoil(klon_glo,nsoilmx,nbsrf)
-cIM "slab" ocean
-      REAL tslab(klon_glo), seaice(klon_glo)
       REAL qsurf(klon_glo,nbsrf)
       REAL qsol(klon_glo)
@@ -492,60 +486,4 @@
       ENDDO
       ENDDO
-c
-cIM "slab" ocean 
-c
-c Lecture de tslab (pour slab ocean seulement):      
-c
-      IF (ocean_in .eq. 'slab  ') then
-        ierr = NF_INQ_VARID (nid, "TSLAB", nvarid)
-        IF (ierr.NE.NF_NOERR) THEN
-          PRINT*, "phyetat0: Le champ <TSLAB> est absent"
-          CALL abort
-        ENDIF
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tslab)
-#else
-        ierr = NF_GET_VAR_REAL(nid, nvarid, tslab)
-#endif
-        IF (ierr.NE.NF_NOERR) THEN
-          PRINT*, "phyetat0: Lecture echouee pour <TSLAB>"
-          CALL abort
-        ENDIF
-        xmin = 1.0E+20
-        xmax = -1.0E+20
-        DO i = 1, klon_glo
-          xmin = MIN(tslab(i),xmin)
-          xmax = MAX(tslab(i),xmax)
-        ENDDO
-        PRINT*,'Min, Max tslab (utilise si OCEAN=slab )', xmin, xmax
-c
-c Lecture de seaice (pour slab ocean seulement):
-c
-        ierr = NF_INQ_VARID (nid, "SEAICE", nvarid)
-        IF (ierr.NE.NF_NOERR) THEN
-          PRINT*, "phyetat0: Le champ <SEAICE> est absent"
-          CALL abort
-        ENDIF
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, seaice)
-#else
-        ierr = NF_GET_VAR_REAL(nid, nvarid, seaice)
-#endif
-        IF (ierr.NE.NF_NOERR) THEN
-          PRINT*, "phyetat0: Lecture echouee pour <SEAICE>"
-          CALL abort
-        ENDIF
-        xmin = 1.0E+20
-        xmax = -1.0E+20
-        DO i = 1, klon_glo
-          xmin = MIN(seaice(i),xmin)
-          xmax = MAX(seaice(i),xmax)
-        ENDDO
-        PRINT*,'Masse de la glace de mer (utilise si OCEAN=slab)', 
-     $  xmin, xmax
-      ELSE
-        tslab = 0.
-        seaice = 0.
-      ENDIF
 c
 c Lecture de l'humidite de l'air juste au dessus du sol:
@@ -1762,6 +1700,4 @@
       call Scatter( wake_fip_glo, wake_fip)
       call Scatter( tsoil,tsoil_p)
-      call Scatter( tslab,tslab_p)
-      call Scatter( seaice,seaice_p)
       call Scatter( qsurf,qsurf_p)
       call Scatter( qsol,qsol_p)
@@ -1796,9 +1732,4 @@
 
 c
-c Initilalize variables in module surface_data
-c
-      ok_veget = ok_veget_in
-      ocean    = ocean_in
-c
 c Initialize module pbl_surface_mod 
 c
@@ -1806,14 +1737,7 @@
      $     evap_p, frugs_p, agesno_p, tsoil_p)
 
-c Initialize ocean module according to ocean type
-      IF ( ocean == 'slab' ) THEN
-c        initilalize module ocean_slab_init
-         CALL ocean_slab_init(dtime, tslab_p, seaice_p, pctsrf)
-      ELSEIF ( ocean == 'couple' ) THEN
-c        initilalize module ocean_cpl_init
+c Initialize module ocean_cpl_mod for the case of coupled ocean
+      IF ( type_ocean == 'couple' ) THEN
          CALL ocean_cpl_init(dtime, rlon, rlat)
-      ELSE
-c        initilalize module ocean_forced_init
-         CALL ocean_forced_init
       ENDIF
 c
Index: LMDZ4/trunk/libf/phylmd/phyredem.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/phyredem.F	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/phyredem.F	(revision 996)
@@ -8,8 +8,6 @@
       USE mod_grid_phy_lmdz
       USE mod_phys_lmdz_para
-      USE ocean_slab_mod,   ONLY : ocean_slab_final
       USE fonte_neige_mod,  ONLY : fonte_neige_final
       USE pbl_surface_mod,  ONLY : pbl_surface_final
-      USE surface_data,     ONLY : ocean, ok_veget
       USE phys_state_var_mod
 
@@ -60,7 +58,5 @@
       REAL wake_fip_glo(klon_glo)
 
-cIM "slab" ocean
       REAL tsoil_p(klon,nsoilmx,nbsrf)
-      REAL tslab_p(klon), seaice_p(klon)
       REAL qsurf_p(klon,nbsrf)
       REAL qsol_p(klon)
@@ -73,5 +69,4 @@
       
       REAL tsoil(klon_glo,nsoilmx,nbsrf)
-      REAL tslab(klon_glo), seaice(klon_glo)
       REAL qsurf(klon_glo,nbsrf)
       REAL qsol(klon_glo)
@@ -103,12 +98,4 @@
 c Get a variable calculated in module fonte_neige_mod
       CALL fonte_neige_final(run_off_lic_0_p)
-
-c If slab ocean then get 2 varaibles from module ocean_slab_mod
-      IF ( ocean == 'slab' ) THEN
-         CALL ocean_slab_final(tslab_p, seaice_p)
-      ELSE 
-         tslab_p(:)  = 0.0
-         seaice_p(:) = 0.0
-      ENDIF      
 
 c======================================================================
@@ -150,6 +137,4 @@
 
       call Gather( tsoil_p,tsoil)
-      call Gather( tslab_p,tslab)
-      call Gather( seaice_p,seaice)
       call Gather( qsurf_p,qsurf)
       call Gather( qsol_p,qsol)
@@ -390,35 +375,4 @@
       ENDDO
 c
-cIM "slab" ocean
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "TSLAB", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "TSLAB", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 33,
-     .                        "Ecart de la SST (pour slab-ocean)")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tslab)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,tslab)
-#endif
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "SEAICE", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "SEAICE", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 33,
-     .                        "Glace de mer kg/m2 (pour slab-ocean)")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,seaice)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,seaice)
-#endif
-c
       DO nsrf = 1, nbsrf
         IF (nsrf.LE.99) THEN
Index: LMDZ4/trunk/libf/phylmd/phys_output_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/phys_output_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/phys_output_mod.F90	(revision 996)
@@ -188,12 +188,5 @@
   integer, dimension(nfiles) , save :: flag_philevsSTD   = (/ 1, 1, 1, 10 /)
 
-  integer, dimension(nfiles) , save :: flag_fluxo        = (/ 1, 1, 10, 10 /)
-  integer, dimension(nfiles) , save :: flag_fluxg        = (/ 1, 1, 10, 10 /)
   integer, dimension(nfiles) , save :: flag_t_oce_sic    = (/ 1, 10, 10, 10 /)
-
-  integer, dimension(nfiles) , save :: flag_lmt_bils     = (/ 1, 1, 10, 10 /)
-  integer, dimension(nfiles) , save :: flag_tslab        = (/ 1, 1, 10, 10 /)
-  integer, dimension(nfiles) , save :: flag_seaice       = (/ 1, 1, 10, 10 /)
-  integer, dimension(nfiles) , save :: flag_siceh        = (/ 1, 1, 10, 10 /)
 
   integer, dimension(nfiles) , save :: flag_weakinv      = (/ 10, 1, 10, 10 /)
@@ -590,18 +583,8 @@
        ENDDO
 
-!IM diagnostiques flux ocean-atm ou ocean-glace de mer
-!IM pour utilisation dans un modele de "slab" ocean
- CALL histdef2d(iff,flag_fluxo,"fluxo","Flux turbulents ocean-atmosphere","W/m2")
- CALL histdef2d(iff,flag_fluxg,"fluxg","Flux turbulents ocean-glace de mer","W/m2")
  CALL histdef2d(iff,flag_t_oce_sic,"t_oce_sic","Temp mixte oce-sic","K")
 
-     IF (OCEAN.EQ.'force ') THEN
- CALL histdef2d(iff,flag_lmt_bils,"lmt_bils","Bilan au sol atmosphere forcee","W/m2")
-     ELSE IF (OCEAN.EQ.'slab  ') THEN
- CALL histdef2d(iff,flag_slab_bils,"slab_bils","Bilan au sol Slab","W/m2") 
- CALL histdef2d(iff,flag_tslab,"tslab", "Slab SST ", "K")
- CALL histdef2d(iff,flag_seaice,"seaice","Slab seaice","kg/m2") 
- CALL histdef2d(iff,flag_siceh,"siceh","Slab seaice height","m")
-     ENDIF
+ IF (ocean=='slab') & 
+      CALL histdef2d(iff,flag_slab_bils, "slab_wbils_oce","Bilan au sol sur ocean slab", "W/m2")
 
  CALL histdef2d(iff,flag_ale_bl,"ale_bl","ALE BL","m2/s2") 
Index: LMDZ4/trunk/libf/phylmd/phys_output_write.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/phys_output_write.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/phys_output_write.h	(revision 996)
@@ -41,5 +41,5 @@
        IF (flag_contfracOR(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),"contfracOR",itau_w,
-     $                   pctsrf_new(:,is_ter))
+     $                   pctsrf(:,is_ter))
        ENDIF
 
@@ -599,29 +599,4 @@
        ENDIF !(bb2.EQ."850".OR.bb2.EQ."700".OR.
        ENDDO
-
-!IM diagnostiques flux ocean-atm ou ocean-glace de mer
-!IM pour utilisation dans un modele de "slab" ocean
-
-      IF (flag_fluxo(iff)<=lev_files(iff)) THEN
-      DO i=1, klon
-       IF (pctsrf(i,is_oce).GT.epsfra) THEN
-        zx_tmp_fi2d(i) = fluxo(i)
-       ELSE
-        zx_tmp_fi2d(i) = 0.
-       ENDIF
-      ENDDO
-      CALL histwrite_phy(nid_files(iff),"fluxo",itau_w,zx_tmp_fi2d)
-      ENDIF
-
-      IF (flag_fluxg(iff)<=lev_files(iff)) THEN
-      DO i=1, klon
-       IF (pctsrf(i,is_sic).GT.epsfra) THEN
-        zx_tmp_fi2d(i) = fluxg(i)
-       ELSE
-        zx_tmp_fi2d(i) = 0.
-       ENDIF
-      ENDDO
-      CALL histwrite_phy(nid_files(iff),"fluxg",itau_w,zx_tmp_fi2d)
-      ENDIF
 
       IF (flag_t_oce_sic(iff)<=lev_files(iff)) THEN
@@ -639,17 +614,5 @@
       ENDIF
 
-      IF (OCEAN.EQ.'force ') THEN
-      IF (flag_lmt_bils(iff)<=lev_files(iff)) THEN
-      DO i=1, klon
-      IF((pctsrf(i,is_oce).GT.epsfra).OR.
-     $   (pctsrf(i,is_sic).GT.epsfra)) THEN
-       zx_tmp_fi2d(i) = (radsol(i) + fluxo(i))*pctsrf(i,is_oce)+
-     $                  fluxg(i)*pctsrf(i,is_sic)
-      ELSE
-       zx_tmp_fi2d(i) = 1.E+20
-      ENDIF
-      ENDDO
-      CALL histwrite_phy(nid_files(iff),"lmt_bils",itau_w,zx_tmp_fi2d)
-      ENDIF
+      IF (type_ocean=='force ') THEN
 
       IF (iflag_coupl.EQ.1) THEN
@@ -683,40 +646,12 @@
        ENDIF
       ENDIF
-
-      ELSE IF (OCEAN.EQ.'slab  ') THEN
-
-      IF (flag_slab_bils(iff)<=lev_files(iff)) THEN
-      DO i=1, klon
-      IF((pctsrf(i,is_oce).GT.epsfra).OR.
-     $   (pctsrf(i,is_sic).GT.epsfra)) THEN
-       zx_tmp_fi2d(i) = (radsol(i) + fluxo(i))*pctsrf(i,is_oce)+
-     $                  fluxg(i)*pctsrf(i,is_sic)
-      ELSE
-       zx_tmp_fi2d(i) = 1.E+20
-      ENDIF
-      ENDDO
-      CALL histwrite_phy(nid_files(iff),"slab_bils",itau_w,zx_tmp_fi2d)
-      ENDIF
-
-      IF (flag_tslab(iff)<=lev_files(iff)) THEN
-      DO i=1, klon
-       IF(pctsrf(i,is_oce).GT.epsfra.OR.
-     $    pctsrf(i,is_sic).GT.epsfra) THEN
-        zx_tmp_fi2d(i)=tslab(i)
-       ELSE
-        zx_tmp_fi2d(i) = 1.E+20
-       ENDIF
-      ENDDO !i=1, klon
-      CALL histwrite_phy(nid_files(iff),"tslab",itau_w,zx_tmp_fi2d)
-      ENDIF
-
-      IF (flag_seaice(iff)<=lev_files(iff)) THEN
-      CALL histwrite_phy(nid_files(iff),"seaice",itau_w,seaice)
-      ENDIF
-
-      IF (flag_siceh(iff)<=lev_files(iff)) THEN
-      CALL histwrite_phy(nid_files(iff),"siceh",itau_w, seaice/1000.)
-      ENDIF
-      ENDIF !OCEAN.EQ.force/slab
+ 
+      ELSE IF (type_ocean=='slab  ') THEN
+
+      IF ( flag_slab_bils(iff)<=lev_files(iff)) 
+     $     CALL histwrite_phy(
+     $     nid_files(iff),"slab_wbils_oce",itau_w,slab_wfbils)
+
+      ENDIF !type_ocean == force/slab
 
       IF (flag_weakinv(iff)<=lev_files(iff)) THEN
Index: LMDZ4/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/physiq.F	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/physiq.F	(revision 996)
@@ -24,12 +24,8 @@
       USE vampir
       USE pbl_surface_mod, ONLY : pbl_surface
-      USE surface_data,     ONLY : ocean, ok_veget
+      USE change_srf_frac_mod
+      USE surface_data,     ONLY : type_ocean, ok_veget
       USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
       USE phys_state_var_mod ! Variables sauvegardees de la physique
-
-
-      USE ocean_slab_mod, ONLY   : ocean_slab_get_vars
-      USE ocean_cpl_mod, ONLY    : ocean_cpl_get_vars
-      USE ocean_forced_mod, ONLY : ocean_forced_get_vars
       USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
       USE phys_output_mod
@@ -120,9 +116,4 @@
       LOGICAL, SAVE :: rnpb=.TRUE.
 c$OMP THREADPRIVATE(rnpb)
-cIM "slab" ocean
-      REAL tslab(klon)    !Temperature du slab-ocean
-      REAL seaice(klon)   !glace de mer (kg/m2) 
-      REAL fluxo(klon)    !flux turbulents ocean-glace de mer 
-      REAL fluxg(klon)    !flux turbulents ocean-atmosphere
       REAL amn, amx
       INTEGER igout
@@ -699,4 +690,6 @@
       REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
 C                             ! type de sous-surface et pondere par la fraction
+      REAL slab_wfbils(klon)  ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
+
       REAL fder(klon)         
       REAL ve(klon) ! integr. verticale du transport meri. de l'energie
@@ -715,8 +708,4 @@
       SAVE lmt_pas                ! frequence de mise a jour
 c$OMP THREADPRIVATE(lmt_pas) 
-cIM
-      REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
-
-cym      SAVE pctsrf                 ! sous-fraction du sol
 
 cIM sorties
@@ -1239,5 +1228,5 @@
 c appel a la lecture du run.def physique
 c
-         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
+         call conf_phys(ok_journe, ok_mensuel,
      .                  ok_instan, ok_hf,
      .                  solarlong0,seuil_inversion,
@@ -1279,5 +1268,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-         CALL phyetat0 ("startphy.nc",ocean, ok_veget,clesphy0,tabcntr0)
+         CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
 cIM begin
           print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
@@ -1288,16 +1277,4 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-       DO i=1,klon
-         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) + 
-     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA) 
-     $       THEN 
-            WRITE(*,*)
-     $   'physiq apres lecture de restart: pb sous surface au point ',
-     $   i, pctsrf(i, 1 : nbsrf)
-         ENDIF 
-      ENDDO
-
 c
 C on remet le calendrier a zero
@@ -1518,5 +1495,5 @@
        call phys_output_open(jjmp1,nqmax,nlevSTD,clevSTD,nbteta, 
      &                        ctetaSTD,dtime,presnivs,ok_veget,
-     &                        ocean,iflag_pbl,ok_mensuel,ok_journe,
+     &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
      &                        ok_hf,ok_instan,nid_files)
 c$OMP END MASTER
@@ -1581,7 +1558,22 @@
 c
       ENDIF
-c
-c   ****************     Fin  de   IF ( debut  )   ***************
-c
+!
+!   ****************     Fin  de   IF ( debut  )   ***************
+!
+!
+! Incrementer le compteur de la physique
+!
+      itap   = itap + 1
+      julien = MOD(NINT(xjour),360)
+      if (julien .eq. 0) julien = 360
+
+!
+! Update fraction of the sub-surfaces (pctsrf) and 
+! initialize, where a new fraction has appeared, all variables depending 
+! on the surface fraction.
+!
+      CALL change_srf_frac(itap, dtime, julien, 
+     *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
+
 
 ! Tendances bidons pour les processus qui n'affectent pas certaines
@@ -1725,10 +1717,4 @@
 cIM END
 c
-c Incrementer le compteur de la physique
-c
-      itap   = itap + 1
-      julien = MOD(NINT(xjour),360)
-      if (julien .eq. 0) julien = 360
-c
 c Mettre en action les conditions aux limites (albedo, sst, etc.).
 c Prescrire l'ozone et calculer l'albedo sur l'ocean.
@@ -1836,5 +1822,5 @@
      s     zxtsol,    zxfluxlat, zt2m,    qsat2m, 
      s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
-     s     ycoefh,    pctsrf_new,                
+     s     ycoefh,    slab_wfbils,                
      d     qsol,      zq2m,      s_pblh,  s_lcl,
      d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
@@ -1847,7 +1833,5 @@
      -     dsens,     devap,     zxsnow,
      -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
-c
-c
-      pctsrf(:,:) = pctsrf_new(:,:)
+
       
 !-----------------------------------------------------------------------------------------
@@ -3332,15 +3316,4 @@
      .     zxfqcalving, zxfqfonte, zxffonte)
 
-      IF (ocean == 'slab') THEN
-         ! Get some variables from module ocean_slab_mod
-         CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
-      ELSEIF (ocean == 'couple') THEN
-         ! Get some variables from module ocean_cpl_mod
-         CALL ocean_cpl_get_vars(fluxo, fluxg)
-      ELSE 
-         ! Get some variables from module ocean_forced_mod
-         CALL ocean_forced_get_vars(fluxo, fluxg)
-      ENDIF
-
 
 c Commente par abderrahmane le 11 2 08
Index: LMDZ4/trunk/libf/phylmd/soil.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/soil.F	(revision 993)
+++ 	(revision )
@@ -1,269 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil,
-     s          pcapcal, pfluxgrd)
-      USE dimphy
-      USE mod_phys_lmdz_para
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteur:  Frederic Hourdin     30/01/92
-c   -------
-c
-c   objet:  computation of : the soil temperature evolution
-c   ------                   the surfacic heat capacity "Capcal"
-c                            the surface conduction flux pcapcal
-c
-c
-c   Method: implicit time integration
-c   -------
-c   Consecutive ground temperatures are related by:
-c           T(k+1) = C(k) + D(k)*T(k)  (1)
-c   the coefficients C and D are computed at the t-dt time-step.
-c   Routine structure:
-c   1)new temperatures are computed  using (1)
-c   2)C and D coefficients are computed from the new temperature
-c     profile for the t+dt time-step
-c   3)the coefficients A and B are computed where the diffusive
-c     fluxes at the t+dt time-step is given by
-c            Fdiff = A + B Ts(t+dt)
-c     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
-c            with F0 = A + B (Ts(t))
-c                 Capcal = B*dt
-c           
-c   Interface:
-c   ----------
-c
-c   Arguments:
-c   ----------
-c   ptimestep            physical timestep (s)
-c   indice               sub-surface index
-c   snow(klon,nbsrf)     snow
-c   ptsrf(klon)          surface temperature at time-step t (K)
-c   ptsoil(klon,nsoilmx) temperature inside the ground (K)
-c   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
-c   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
-c   
-c=======================================================================
-c   declarations:
-c   -------------
-
-cym#include "dimensions.h"
-#include "YOMCST.h"
-cym#include "dimphy.h"
-#include "dimsoil.h"
-#include "indicesol.h"
-#include "comsoil.h"
-
-c-----------------------------------------------------------------------
-c  arguments
-c  ---------
-
-      REAL ptimestep
-      INTEGER indice, knon
-      REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)
-      REAL pcapcal(klon),pfluxgrd(klon)
-
-c-----------------------------------------------------------------------
-c  local arrays
-c  ------------
-
-      INTEGER ig,jk
-c@$$      REAL zdz2(nsoilmx),z1(klon)
-      REAL zdz2(nsoilmx),z1(klon,nbsrf)
-      REAL,SAVE :: min_period,dalph_soil
-c$OMP THREADPRIVATE( min_period,dalph_soil)     
-      REAL ztherm_i(klon)
-
-c   local saved variables:
-c   ----------------------
-      REAL dz1(nsoilmx),dz2(nsoilmx)
-c@$$          REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
-cym      REAL zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf)
-      REAL,ALLOCATABLE,SAVE ::  zc(:,:,:),zd(:,:,:)
-c$OMP THREADPRIVATE(zc,zd)
-      REAL lambda
-cym      SAVE dz1,dz2,zc,zd,lambda
-      SAVE dz1,dz2,lambda
-c$OMP THREADPRIVATE(dz1,dz2,lambda)
-      LOGICAL firstcall, firstsurf(nbsrf)
-      SAVE firstcall, firstsurf
-c$OMP THREADPRIVATE(firstcall, firstsurf)
-      DATA firstcall/.true./
-      DATA firstsurf/.TRUE.,.TRUE.,.TRUE.,.TRUE./
-
-      LOGICAL,SAVE :: First=.true.
-c$OMP THREADPRIVATE(First)
-c-----------------------------------------------------------------------
-c   Depthts:
-c   --------
-
-      REAL fz,rk,fz1,rk1,rk2
-      fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
-      pfluxgrd(:) = 0.
-c   calcul de l'inertie thermique a partir de la variable rnat.
-c   on initialise a inertie_ice meme au-dessus d'un point de mer au cas 
-c   ou le point de mer devienne point de glace au pas suivant
-c   on corrige si on a un point de terre avec ou sans glace
-c
-      IF (first) THEN
-        allocate(zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf))
-	first=.false.
-      ENDIF
-      
-      IF (indice.EQ.is_sic) THEN
-         DO ig = 1, knon
-            ztherm_i(ig)   = inertie_ice
-            IF (snow(ig).GT.0.0) ztherm_i(ig)   = inertie_sno
-         ENDDO
-      ELSE IF (indice.EQ.is_lic) THEN
-         DO ig = 1, knon
-            ztherm_i(ig)   = inertie_ice
-            IF (snow(ig).GT.0.0) ztherm_i(ig)   = inertie_sno
-         ENDDO
-      ELSE IF (indice.EQ.is_ter) THEN
-         DO ig = 1, knon
-            ztherm_i(ig)   = inertie_sol
-            IF (snow(ig).GT.0.0) ztherm_i(ig)   = inertie_sno
-         ENDDO
-      ELSE IF (indice.EQ.is_oce) THEN
-         DO ig = 1, knon
-            ztherm_i(ig)   = inertie_ice
-         ENDDO
-      ELSE
-         PRINT*, "valeur d indice non prevue", indice
-         CALL abort
-      ENDIF
-
-
-c@$$      IF (firstcall) THEN
-      IF (firstsurf(indice)) THEN 
-
-c-----------------------------------------------------------------------
-c   ground levels 
-c   grnd=z/l where l is the skin depth of the diurnal cycle:
-c   --------------------------------------------------------
-
-         min_period=1800. ! en secondes
-         dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
-c$OMP MASTER
-         IF (is_mpi_root) THEN
-         OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
-         READ(99,*) min_period
-         READ(99,*) dalph_soil
-         PRINT*,'Discretization for the soil model'
-         PRINT*,'First level e-folding depth',min_period,
-     s   '   dalph',dalph_soil
-         CLOSE(99)
-9999     CONTINUE
-         ENDIF  ! is_mpi_root
-c$OMP END MASTER
-         CALL bcast(min_period)
-	 CALL bcast(dalph_soil)
-
-c   la premiere couche represente un dixieme de cycle diurne
-         fz1=sqrt(min_period/3.14)
-
-         DO jk=1,nsoilmx
-            rk1=jk
-            rk2=jk-1
-            dz2(jk)=fz(rk1)-fz(rk2)
-         ENDDO
-         DO jk=1,nsoilmx-1
-            rk1=jk+.5
-            rk2=jk-.5
-            dz1(jk)=1./(fz(rk1)-fz(rk2))
-         ENDDO
-         lambda=fz(.5)*dz1(1)
-         PRINT*,'full layers, intermediate layers (seconds)'
-         DO jk=1,nsoilmx
-            rk=jk
-            rk1=jk+.5
-            rk2=jk-.5
-            PRINT *,'fz=',
-     .               fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
-         ENDDO
-C PB
-         firstsurf(indice) = .FALSE. 
-c@$$         firstcall =.false.
-
-c   Initialisations:
-c   ----------------
-
-      ELSE   !--not firstcall
-c-----------------------------------------------------------------------
-c   Computation of the soil temperatures using the Cgrd and Dgrd
-c  coefficient computed at the previous time-step:
-c  -----------------------------------------------
-
-c    surface temperature
-         DO ig=1,knon
-            ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/
-     s      (lambda*(1.-zd(ig,1,indice))+1.)
-         ENDDO
-
-c   other temperatures
-         DO jk=1,nsoilmx-1
-            DO ig=1,knon
-               ptsoil(ig,jk+1)=zc(ig,jk,indice)+zd(ig,jk,indice)
-     $            *ptsoil(ig,jk)
-            ENDDO
-         ENDDO
-
-      ENDIF !--not firstcall
-c-----------------------------------------------------------------------
-c   Computation of the Cgrd and Dgrd coefficient for the next step:
-c   ---------------------------------------------------------------
-
-c@$$  PB ajout pour cas glace de mer
-      IF (indice .EQ. is_sic) THEN
-          DO ig = 1 , knon
-            ptsoil(ig,nsoilmx) = RTT - 1.8
-          END DO 
-      ENDIF 
-
-      DO jk=1,nsoilmx
-         zdz2(jk)=dz2(jk)/ptimestep
-      ENDDO
-
-      DO ig=1,knon
-         z1(ig,indice)=zdz2(nsoilmx)+dz1(nsoilmx-1)
-         zc(ig,nsoilmx-1,indice)=
-     $       zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig,indice)
-         zd(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1(ig,indice)
-      ENDDO
-
-      DO jk=nsoilmx-1,2,-1
-         DO ig=1,knon
-            z1(ig,indice)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)
-     $         *(1.-zd(ig,jk,indice)))
-            zc(ig,jk-1,indice)=
-     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice))
-     $          *z1(ig,indice)
-            zd(ig,jk-1,indice)=dz1(jk-1)*z1(ig,indice)
-         ENDDO
-      ENDDO
-
-c-----------------------------------------------------------------------
-c   computation of the surface diffusive flux from ground and
-c   calorific capacity of the ground:
-c   ---------------------------------
-
-      DO ig=1,knon
-         pfluxgrd(ig)=ztherm_i(ig)*dz1(1)*
-     s   (zc(ig,1,indice)+(zd(ig,1,indice)-1.)*ptsoil(ig,1))
-         pcapcal(ig)=ztherm_i(ig)*
-     s   (dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))
-         z1(ig,indice)=lambda*(1.-zd(ig,1,indice))+1.
-         pcapcal(ig)=pcapcal(ig)/z1(ig,indice)
-         pfluxgrd(ig) = pfluxgrd(ig)
-     s   + pcapcal(ig) * (ptsoil(ig,1) * z1(ig,indice)
-     $       - lambda * zc(ig,1,indice)
-     $       - ptsrf(ig))
-     s   /ptimestep
-      ENDDO
-
-      RETURN
-      END
Index: LMDZ4/trunk/libf/phylmd/soil.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/soil.F90	(revision 996)
+++ LMDZ4/trunk/libf/phylmd/soil.F90	(revision 996)
@@ -0,0 +1,278 @@
+!
+! $Header$
+!
+SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, &
+     ptsoil, pcapcal, pfluxgrd)
+  
+  USE dimphy
+  USE mod_phys_lmdz_para
+  IMPLICIT NONE
+
+!=======================================================================
+!
+!   Auteur:  Frederic Hourdin     30/01/92
+!   -------
+!
+!   Object:  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)  (*)
+!   The coefficients C and D are computed at the t-dt time-step.
+!   Routine structure:
+!   1) C and D coefficients are computed from the old temperature
+!   2) new temperatures are computed using (*)
+!   3) C and D coefficients are computed from the new temperature
+!      profile for the t+dt time-step
+!   4) 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:
+!   ----------
+!   ptimestep            physical timestep (s)
+!   indice               sub-surface index
+!   snow(klon)           snow
+!   ptsrf(klon)          surface temperature at time-step t (K)
+!   ptsoil(klon,nsoilmx) temperature inside the ground (K)
+!   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
+!   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
+!   
+!=======================================================================
+  INCLUDE "YOMCST.h"
+  INCLUDE "dimsoil.h"
+  INCLUDE "indicesol.h"
+  INCLUDE "comsoil.h"
+!-----------------------------------------------------------------------
+! Arguments
+! ---------
+  REAL, INTENT(IN)                     :: ptimestep
+  INTEGER, INTENT(IN)                  :: indice, knon
+  REAL, DIMENSION(klon), INTENT(IN)    :: snow
+  REAL, DIMENSION(klon), INTENT(IN)    :: ptsrf
+  
+  REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil
+  REAL, DIMENSION(klon), INTENT(OUT)           :: pcapcal
+  REAL, DIMENSION(klon), INTENT(OUT)           :: pfluxgrd
+
+!-----------------------------------------------------------------------
+! Local variables
+! ---------------
+  INTEGER                             :: ig, jk, ierr
+  REAL                                :: min_period,dalph_soil
+  REAL, DIMENSION(nsoilmx)            :: zdz2
+  REAL                                :: z1s
+  REAL, DIMENSION(klon)               :: ztherm_i
+  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef
+
+! Local saved variables
+! ---------------------
+  REAL, SAVE                     :: lambda
+!$OMP THREADPRIVATE(lambda)
+  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
+!$OMP THREADPRIVATE(dz1,dz2)
+  LOGICAL, SAVE                  :: firstcall=.TRUE.
+!$OMP THREADPRIVATE(firstcall)
+    
+!-----------------------------------------------------------------------
+!   Depthts:
+!   --------
+  REAL fz,rk,fz1,rk1,rk2
+  fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
+
+
+!-----------------------------------------------------------------------
+! Calculation of some constants
+! NB! These constants do not depend on the sub-surfaces
+!-----------------------------------------------------------------------
+
+  IF (firstcall) THEN
+!-----------------------------------------------------------------------
+!   ground levels 
+!   grnd=z/l where l is the skin depth of the diurnal cycle:
+!-----------------------------------------------------------------------
+
+     min_period=1800. ! en secondes
+     dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
+!$OMP MASTER
+     IF (is_mpi_root) THEN
+        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
+        IF (ierr == 0) THEN ! Read file only if it exists
+           READ(99,*) min_period
+           READ(99,*) dalph_soil
+           PRINT*,'Discretization for the soil model'
+           PRINT*,'First level e-folding depth',min_period, &
+                '   dalph',dalph_soil
+           CLOSE(99)
+        END IF
+     ENDIF
+!$OMP END MASTER
+     CALL bcast(min_period)
+     CALL bcast(dalph_soil)
+
+!   la premiere couche represente un dixieme de cycle diurne
+     fz1=SQRT(min_period/3.14)
+     
+     DO jk=1,nsoilmx
+        rk1=jk
+        rk2=jk-1
+        dz2(jk)=fz(rk1)-fz(rk2)
+     ENDDO
+     DO jk=1,nsoilmx-1
+        rk1=jk+.5
+        rk2=jk-.5
+        dz1(jk)=1./(fz(rk1)-fz(rk2))
+     ENDDO
+     lambda=fz(.5)*dz1(1)
+     PRINT*,'full layers, intermediate layers (seconds)'
+     DO jk=1,nsoilmx
+        rk=jk
+        rk1=jk+.5
+        rk2=jk-.5
+        PRINT *,'fz=', &
+             fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
+     ENDDO
+
+     firstcall =.FALSE.
+  END IF
+
+
+!-----------------------------------------------------------------------
+!   Calcul de l'inertie thermique a partir de la variable rnat.
+!   on initialise a inertie_ice meme au-dessus d'un point de mer au cas 
+!   ou le point de mer devienne point de glace au pas suivant
+!   on corrige si on a un point de terre avec ou sans glace
+!
+!-----------------------------------------------------------------------
+  IF (indice == is_sic) THEN
+     DO ig = 1, knon
+        ztherm_i(ig)   = inertie_ice
+        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
+     ENDDO
+  ELSE IF (indice == is_lic) THEN
+     DO ig = 1, knon
+        ztherm_i(ig)   = inertie_ice
+        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
+     ENDDO
+  ELSE IF (indice == is_ter) THEN
+     DO ig = 1, knon
+        ztherm_i(ig)   = inertie_sol
+        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
+     ENDDO
+  ELSE IF (indice == is_oce) THEN
+     DO ig = 1, knon
+        ztherm_i(ig)   = inertie_ice
+     ENDDO
+  ELSE
+     PRINT*, "valeur d indice non prevue", indice
+     CALL abort
+  ENDIF
+
+
+!-----------------------------------------------------------------------
+! 1)
+! Calculation of Cgrf and Dgrd coefficients using soil temperature from 
+! previous time step.
+!
+! These variables are recalculated on the local compressed grid instead 
+! of saved in restart file.
+!-----------------------------------------------------------------------
+  DO jk=1,nsoilmx
+     zdz2(jk)=dz2(jk)/ptimestep
+  ENDDO
+  
+  DO ig=1,knon
+     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
+     C_coef(ig,nsoilmx-1,indice)= &
+          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
+     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
+  ENDDO
+  
+  DO jk=nsoilmx-1,2,-1
+     DO ig=1,knon
+        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
+             *(1.-D_coef(ig,jk,indice)))
+        C_coef(ig,jk-1,indice)= &
+             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
+        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
+     ENDDO
+  ENDDO
+
+!-----------------------------------------------------------------------
+! 2)
+! Computation of the soil temperatures using the Cgrd and Dgrd
+! coefficient computed above
+!
+!-----------------------------------------------------------------------
+
+!    Surface temperature
+  DO ig=1,knon
+     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
+          (lambda*(1.-D_coef(ig,1,indice))+1.)
+  ENDDO
+  
+!   Other temperatures
+  DO jk=1,nsoilmx-1
+     DO ig=1,knon
+        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
+             *ptsoil(ig,jk)
+     ENDDO
+  ENDDO
+
+  IF (indice == is_sic) THEN
+     DO ig = 1 , knon
+        ptsoil(ig,nsoilmx) = RTT - 1.8
+     END DO
+  ENDIF
+
+!-----------------------------------------------------------------------
+! 3)
+! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil 
+! temperature
+!-----------------------------------------------------------------------
+  DO ig=1,knon
+     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
+     C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
+     D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
+  ENDDO
+  
+  DO jk=nsoilmx-1,2,-1
+     DO ig=1,knon
+        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
+             *(1.-D_coef(ig,jk,indice)))
+        C_coef(ig,jk-1,indice) = &
+             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
+        D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s
+     ENDDO
+  ENDDO
+
+!-----------------------------------------------------------------------
+! 4)
+! Computation of the surface diffusive flux from ground and
+! calorific capacity of the ground
+!-----------------------------------------------------------------------
+  DO ig=1,knon
+     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
+          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
+     pcapcal(ig)  = ztherm_i(ig)* &
+          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
+     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
+     pcapcal(ig)  = pcapcal(ig)/z1s
+     pfluxgrd(ig) = pfluxgrd(ig) &
+          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
+          - lambda * C_coef(ig,1,indice) &
+          - ptsrf(ig)) &
+          /ptimestep
+  ENDDO
+    
+END SUBROUTINE soil
Index: LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90	(revision 996)
@@ -8,10 +8,4 @@
 ! This module is used when no external land model is choosen.
 !
-  USE fonte_neige_mod
-  USE calcul_fluxs_mod
-  USE dimphy
-  USE mod_grid_phy_lmdz
-  USE mod_phys_lmdz_para
- 
   IMPLICIT NONE
 
@@ -26,4 +20,12 @@
        fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l)
 
+    USE limit_read_mod
+    USE surface_data
+    USE fonte_neige_mod
+    USE calcul_fluxs_mod
+    USE cpl_mod
+    USE dimphy
+    USE mod_grid_phy_lmdz
+    USE mod_phys_lmdz_para
 !****************************************************************************************
 ! Bucket calculations for surface. 
@@ -74,4 +76,5 @@
     REAL, DIMENSION(klon) :: zfra
     REAL, DIMENSION(klon) :: radsol       ! total net radiance at surface
+    REAL, DIMENSION(klon) :: dummy_riverflow,dummy_coastalflow 
     INTEGER               :: i
 !
@@ -82,8 +85,7 @@
 !* Read from limit.nc : albedo(alb_lim), length of rugosity(z0_new)
 !
-    CALL interfsur_lim(itime, dtime, jour, &
-         knon, knindex, debut,  &
-         alb_lim, z0_new)
-    
+    CALL limit_read_rug_alb(itime, dtime, jour,&
+         knon, knindex, &
+         z0_new, alb_lim)
 !
 !* Calcultaion of fluxes 
@@ -145,6 +147,18 @@
 !* Calculate the rugosity
 !
-    z0_new = SQRT(z0_new**2+rugoro**2)
-       
+    DO i = 1, knon
+       z0_new(i) = SQRT(z0_new(i)**2+rugoro(i)**2)
+    END DO
+
+!* Send to coupler
+!  The run-off from river and coast are not calculated in the bucket modele.
+!  For testing purpose of the coupled modele we put the run-off to zero.
+    IF (type_ocean=='couple') THEN
+       dummy_riverflow(:)   = 0.0
+       dummy_coastalflow(:) = 0.0
+       CALL cpl_send_land_fields(itime, knon, knindex, &
+            dummy_riverflow, dummy_coastalflow)
+    ENDIF
+
 !
 !* End
@@ -154,176 +168,3 @@
 !****************************************************************************************
 !
-  SUBROUTINE interfsur_lim(itime, dtime, jour, &
-       knon, knindex, debut, &
-       lmt_alb_p, lmt_rug_p)
-    
-! Cette routine sert d'interface entre le modele atmospherique et un fichier
-! de conditions aux limites
-!
-! L. Fairhead 02/2000
-!
-! input:
-!   itime        numero du pas de temps courant
-!   dtime        pas de temps de la physique (en s)
-!   jour         jour a lire dans l'annee
-!   knon         nombre de points dans le domaine a traiter
-!   knindex      index des points de la surface a traiter
-!   debut        logical: 1er appel a la physique (initialisation)
-!
-! output:
-!   lmt_alb_p      Albedo lu 
-!   lmt_rug_p      longueur de rugosite lue
-
-    INCLUDE "netcdf.inc"
-
-! Input variables
-!****************************************************************************************
-    INTEGER, INTENT(IN)                      :: itime
-    REAL   , INTENT(IN)                      :: dtime
-    INTEGER, INTENT(IN)                      :: jour
-    INTEGER, INTENT(IN)                      :: knon
-    INTEGER, DIMENSION(klon_loc), INTENT(IN) :: knindex
-    LOGICAL, INTENT(IN)                      :: debut
-
-! Output variables
-!****************************************************************************************
-    REAL, INTENT(out), DIMENSION(klon_loc)   :: lmt_alb_p
-    REAL, INTENT(out), DIMENSION(klon_loc)   :: lmt_rug_p
-
-! Local variables with attribute SAVE
-!****************************************************************************************
-    INTEGER,SAVE   :: lmt_pas     ! frequence de lecture des conditions limites 
-                                  ! (en pas de physique)
-    !$OMP THREADPRIVATE(lmt_pas) 
-    LOGICAL,SAVE   :: deja_lu_sur ! pour indiquer que le jour a lire a deja
-                                  ! lu pour une surface precedente
-    !$OMP THREADPRIVATE(deja_lu_sur) 
-    INTEGER,SAVE                           :: jour_lu_sur 
-    !$OMP THREADPRIVATE(jour_lu_sur) 
-    CHARACTER (len = 20),SAVE              :: fich ='limit.nc'
-    !$OMP THREADPRIVATE(fich) 
-    LOGICAL,SAVE                           :: check = .FALSE.
-    !$OMP THREADPRIVATE(check) 
-! Champs lus dans le fichier de CL
-    REAL, ALLOCATABLE , SAVE, DIMENSION(:) :: alb_lu_p, rug_lu_p
-    !$OMP THREADPRIVATE(alb_lu_p, rug_lu_p) 
-
-! quelques variables pour netcdf
-    INTEGER ,SAVE                          :: nid, nvarid
-    !$OMP THREADPRIVATE(nid, nvarid) 
-    INTEGER, DIMENSION(2),SAVE             :: start, epais
-    !$OMP THREADPRIVATE(start, epais) 
-
-! Other local variables
-!****************************************************************************************
-    INTEGER                                :: ii, ierr
-    CHARACTER (len = 20)                   :: modname = 'interfsur_lim'
-    CHARACTER (len = 80)                   :: abort_message
-    REAL, DIMENSION(klon_glo)              :: alb_lu
-    REAL, DIMENSION(klon_glo)              :: rug_lu
-
-!
-! End delaration
-!****************************************************************************************
-
-    IF (debut) THEN
-       lmt_pas = NINT(86400./dtime * 1.0) ! pour une lecture une fois par jour
-       jour_lu_sur = jour - 1
-       ALLOCATE(alb_lu_p(klon_loc))
-       ALLOCATE(rug_lu_p(klon_loc))
-    ENDIF
-    
-    IF ((jour - jour_lu_sur) /= 0) deja_lu_sur = .FALSE.
-  
-    IF (check) WRITE(*,*) modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur 
-    IF (check) WRITE(*,*) modname,':: itime, lmt_pas', itime, lmt_pas
-    IF (check) CALL flush(6)
-
-!    
-! Tester d'abord si c'est le moment de lire le fichier
-!
-    IF (MOD(itime-1, lmt_pas) == 0 .AND. .NOT. deja_lu_sur) THEN
-
-!
-! Ouverture et lecture du fichier
-!
-!$OMP MASTER
-       IF (is_mpi_root) THEN
-          fich = TRIM(fich)
-          IF (check) WRITE(*,*)modname,' ouverture fichier ',fich
-          IF (check) CALL flush(6)
-          ierr = NF_OPEN (fich, NF_NOWRITE,nid)
-          IF (ierr.NE.NF_NOERR) THEN
-             abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
-             CALL abort_gcm(modname,abort_message,1)
-          ENDIF
-! 
-! La tranche de donnees a lire:
-          start(1) = 1
-          start(2) = jour
-          epais(1) = klon_glo
-          epais(2) = 1
-!
-! Lecture albedo
-          ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
-          IF (ierr /= NF_NOERR) THEN
-             abort_message = 'Le champ <ALB> est absent'
-             CALL abort_gcm(modname,abort_message,1)
-          ENDIF
-#ifdef NC_DOUBLE
-          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
-#else
-          ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
-#endif
-          IF (ierr /= NF_NOERR) THEN
-             abort_message = 'Lecture echouee pour <ALB>'
-             CALL abort_gcm(modname,abort_message,1)
-          ENDIF
-!
-! Lecture rugosite!
-          ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
-          IF (ierr /= NF_NOERR) THEN
-             abort_message = 'Le champ <RUG> est absent'
-             CALL abort_gcm(modname,abort_message,1)
-          ENDIF
-#ifdef NC_DOUBLE
-          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
-#else
-          ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
-#endif
-          IF (ierr /= NF_NOERR) THEN
-             abort_message = 'Lecture echouee pour <RUG>'
-             CALL abort_gcm(modname,abort_message,1)
-          ENDIF
-
-!
-! Fin de lecture
-          ierr = NF_CLOSE(nid)
-
-       ENDIF ! is_mpi_root
-!$OMP END MASTER
-
-       CALL Scatter(alb_lu,alb_lu_p)
-       CALL Scatter(rug_lu,rug_lu_p)
-
-       deja_lu_sur = .TRUE.
-       jour_lu_sur = jour       
-
-    ENDIF
-  
-!
-! Recopie des variables dans les champs de sortie
-!
-    lmt_alb_p(:) = 999999.
-    lmt_rug_p(:) = 999999.
-    DO ii = 1, knon
-       lmt_alb_p(ii) = alb_lu_p(knindex(ii))
-       lmt_rug_p(ii) = rug_lu_p(knindex(ii))
-    ENDDO
-    
-
-  END SUBROUTINE interfsur_lim
-!
-!****************************************************************************************
-!
 END MODULE surf_land_bucket_mod
Index: LMDZ4/trunk/libf/phylmd/surf_land_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_land_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/surf_land_mod.F90	(revision 996)
@@ -26,5 +26,5 @@
        snow, qsol, agesno, tsoil, &
        z0_new, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-       qsurf, tsurf_new, dflux_s, dflux_l, pctsrf_ter, &
+       qsurf, tsurf_new, dflux_s, dflux_l, &
        lwdown_m)
 
@@ -75,5 +75,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_ter
 
 ! Local variables
@@ -154,11 +153,4 @@
     ENDIF ! ok_veget
 
-!****************************************************************************************
-! Return the pourcentage of land in each grid cell, even if not changed in here!
-!
-!****************************************************************************************
-    pctsrf_ter(:) = pctsrf(:,is_ter)
-
-
   END SUBROUTINE surf_land
 !
Index: LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 996)
@@ -15,5 +15,5 @@
   USE intersurf     ! module d'ORCHIDEE
   USE cpl_mod,      ONLY : cpl_send_land_fields
-  USE surface_data, ONLY : ocean, ok_veget
+  USE surface_data, ONLY : type_ocean
   USE comgeomphy,   ONLY : cuphy, cvphy
   USE mod_grid_phy_lmdz
@@ -199,5 +199,4 @@
 
     IF (check) WRITE(lunout,*)'Entree ', modname
-    IF (check) WRITE(lunout,*)'ok_veget = ',ok_veget
   
 ! Initialisation
@@ -416,5 +415,5 @@
 !* Send to coupler
 !
-    IF (ocean=='couple') THEN
+    IF (type_ocean=='couple') THEN
        CALL cpl_send_land_fields(itime, knon, knindex, &
             riverflow, coastalflow)
Index: LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90	(revision 996)
@@ -5,5 +5,5 @@
   
   USE dimphy
-  USE surface_data,     ONLY : ocean, calice, calsno
+  USE surface_data,     ONLY : type_ocean, calice, calsno
   USE fonte_neige_mod,  ONLY : fonte_neige, run_off_lic
   USE cpl_mod,          ONLY : cpl_send_landice_fields
@@ -23,5 +23,5 @@
        snow, qsurf, qsol, agesno, &
        tsoil, z0_new, alb1, alb2, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_lic)
+       tsurf_new, dflux_s, dflux_l)
 
     INCLUDE "indicesol.h"
@@ -64,5 +64,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_lic
 
 ! Local variables
@@ -153,12 +152,4 @@
     z0_new(:) = rugoro(:)
 
-
-!****************************************************************************************
-! Return the pourcentage for this sub-surface
-!
-!****************************************************************************************
-    pctsrf_lic(:) = pctsrf(:,is_lic)
-    
-
 !****************************************************************************************
 ! Send run-off on land-ice to coupler if coupled ocean.
@@ -166,5 +157,5 @@
 !
 !****************************************************************************************
-    IF (ocean=='couple') THEN
+    IF (type_ocean=='couple') THEN
        CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic)
     ENDIF
Index: LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 996)
@@ -5,5 +5,5 @@
 
   USE dimphy
-  USE surface_data, ONLY     : ocean
+  USE surface_data, ONLY     : type_ocean
   USE ocean_forced_mod, ONLY : ocean_forced_noice
   USE ocean_slab_mod, ONLY   : ocean_slab_noice
@@ -17,5 +17,5 @@
 !
   SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
-       rugos, windsp, rmu0, fder, &
+       rugos, windsp, rmu0, fder, tsurf_in, &
        itime, dtime, jour, knon, knindex, &
        debut, &
@@ -25,5 +25,5 @@
        snow, qsurf, agesno, &
        z0_new, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+       tsurf_new, dflux_s, dflux_l, lmt_bils)
 !
 ! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force, 
@@ -47,4 +47,5 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: rmu0  
     REAL, DIMENSION(klon), INTENT(IN)        :: fder
+    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
     REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
     REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
@@ -74,6 +75,5 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
-
+    REAL, DIMENSION(klon), INTENT(OUT)       :: lmt_bils
 
 ! Local variables
@@ -97,5 +97,5 @@
 ! Switch according to type of ocean (couple, slab or forced)
 !****************************************************************************************
-    SELECT CASE(ocean)
+    SELECT CASE(type_ocean)
     CASE('couple')
        CALL ocean_cpl_noice( &
@@ -106,18 +106,18 @@
             p1lay, tq_cdrag, precip_rain, precip_snow,temp_air,spechum,&
             petAcoef, peqAcoef, petBcoef, peqBcoef, &
-            ps, u1_lay, v1_lay, pctsrf, &
+            ps, u1_lay, v1_lay, &
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, &
-            tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+            tsurf_new, dflux_s, dflux_l)
 
     CASE('slab')
        CALL ocean_slab_noice( &
-            dtime, knon, knindex, &
+            itime, dtime, jour, knon, knindex, &
             p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
             petAcoef, peqAcoef, petBcoef, peqBcoef, &
-            ps, u1_lay, v1_lay, &
+            ps, u1_lay, v1_lay, tsurf_in, &
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, &
-            tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+            tsurf_new, dflux_s, dflux_l, lmt_bils)
        
     CASE('force')
@@ -131,5 +131,5 @@
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, &
-            tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+            tsurf_new, dflux_s, dflux_l)
     END SELECT
 
Index: LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 996)
@@ -5,6 +5,5 @@
 
   USE dimphy
-  USE surface_data, ONLY     : ocean
-  USE ocean_slab_mod, ONLY   : ocean_slab_ice
+  USE surface_data
   USE ocean_forced_mod, ONLY : ocean_forced_ice
   USE ocean_cpl_mod, ONLY    : ocean_cpl_ice
@@ -24,5 +23,5 @@
        snow, qsurf, qsol, agesno, tsoil, &
        z0_new, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
+       tsurf_new, dflux_s, dflux_l)
 !
 ! This subroutine will make a call to ocean_XXX_ice according to the ocean mode (force, 
@@ -70,5 +69,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_sic
 
 ! Local arguments
@@ -91,39 +89,40 @@
 !
 !****************************************************************************************
-    SELECT CASE(ocean)
-    CASE('couple')
+    IF (type_ocean == 'couple') THEN
+       
        CALL ocean_cpl_ice( &
-          rlon, rlat, swnet, lwnet, alb1, & 
-          fder, & 
-          itime, dtime, knon, knindex, &
-          lafin,&
-          p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
-          petAcoef, peqAcoef, petBcoef, peqBcoef, &
-          ps, u1_lay, v1_lay, pctsrf, &
-          radsol, snow, qsurf, &
-          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-          tsurf_new, dflux_s, dflux_l, pctsrf_sic)
+            rlon, rlat, swnet, lwnet, alb1, & 
+            fder, & 
+            itime, dtime, knon, knindex, &
+            lafin,&
+            p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
+            petAcoef, peqAcoef, petBcoef, peqBcoef, &
+            ps, u1_lay, v1_lay, pctsrf, &
+            radsol, snow, qsurf, &
+            alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+            tsurf_new, dflux_s, dflux_l)
+       
+    ELSE IF (type_ocean == 'force' .OR. (type_ocean == 'slab' .AND. version_ocean=='sicOBS')) THEN
+       CALL ocean_forced_ice(itime, dtime, jour, knon, knindex, &
+            debut, &
+            tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
+            petAcoef, peqAcoef, petBcoef, peqBcoef, &
+            ps, u1_lay, v1_lay, &
+            radsol, snow, qsol, agesno, tsoil, &
+            qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+            tsurf_new, dflux_s, dflux_l)
 
-    CASE('slab')
-       CALL ocean_slab_ice( & 
-          itime, dtime, jour, knon, knindex, &
-          debut, &
-          tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
-          petAcoef, peqAcoef, petBcoef, peqBcoef, &
-          ps, u1_lay, v1_lay, &
-          radsol, snow, qsurf, qsol, agesno, tsoil, &
-          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-          tsurf_new, dflux_s, dflux_l, pctsrf_sic)
-    
-    CASE('force')
-       CALL ocean_forced_ice(itime, dtime, jour, knon, knindex, &
-          debut, &
-          tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
-          petAcoef, peqAcoef, petBcoef, peqBcoef, &
-          ps, u1_lay, v1_lay, &
-          radsol, snow, qsol, agesno, tsoil, &
-          qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-          tsurf_new, dflux_s, dflux_l, pctsrf_sic)
-    END SELECT
+    ELSE IF (type_ocean == 'slab') THEN
+!!$       CALL ocean_slab_ice( & 
+!!$          itime, dtime, jour, knon, knindex, &
+!!$          debut, &
+!!$          tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
+!!$          petAcoef, peqAcoef, petBcoef, peqBcoef, &
+!!$          ps, u1_lay, v1_lay, pctsrf, &
+!!$          radsol, snow, qsurf, qsol, agesno, tsoil, &
+!!$          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+!!$          tsurf_new, dflux_s, dflux_l)
+
+    END IF
 
 !****************************************************************************************
Index: LMDZ4/trunk/libf/phylmd/surface_data.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surface_data.F90	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/surface_data.F90	(revision 996)
@@ -4,14 +4,18 @@
 MODULE surface_data
 
-  USE dimphy, ONLY : klon
-
   REAL, PARAMETER        :: calice=1.0/(5.1444e+06*0.15)
   REAL, PARAMETER        :: tau_gl=86400.*5.
   REAL, PARAMETER        :: calsno=1./(2.3867e+06*.15)
   
-  LOGICAL, SAVE          :: ok_veget       ! true for use of vegetation model ORCHIDEE
+  LOGICAL, SAVE          :: ok_veget      ! true for use of vegetation model ORCHIDEE
   !$OMP THREADPRIVATE(ok_veget)
-  CHARACTER(len=6), SAVE :: ocean          ! force/slab/couple
-  !$OMP THREADPRIVATE(ocean)
+
+  CHARACTER(len=6), SAVE :: type_ocean    ! force/slab/couple
+  !$OMP THREADPRIVATE(type_ocean)
+
+  ! if type_ocean=couple : version_ocean=opa8 ou nemo
+  ! if type_ocean=slab   : version_ocean=sicOBS
+  CHARACTER(len=6), SAVE :: version_ocean 
+  !$OMP THREADPRIVATE(version_ocean)
 
 END MODULE surface_data
Index: LMDZ4/trunk/libf/phylmd/write_histday.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_histday.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/write_histday.h	(revision 996)
@@ -32,7 +32,7 @@
       CALL histwrite_phy(nid_day,"contfracATM",itau_w,zx_tmp_fi2d)
 c
-cym      CALL gr_fi_ecrit(1,klon,iim,jjmp1,pctsrf_new(:,is_ter),zx_tmp_2d)
+cym      CALL gr_fi_ecrit(1,klon,iim,jjmp1,pctsrf(:,is_ter),zx_tmp_2d)
       CALL histwrite_phy(nid_day,"contfracOR",itau_w,
-     &                   pctsrf_new(:,is_ter))
+     &                   pctsrf(:,is_ter))
 c
 
Index: LMDZ4/trunk/libf/phylmd/write_histday_seri.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_histday_seri.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/write_histday_seri.h	(revision 996)
@@ -191,5 +191,5 @@
 c
 c     ok_msk=.TRUE.
-c     msk(1:klon)=pctsrf_new(1:klon,is_ter)
+c     msk(1:klon)=pctsrf(1:klon,is_ter)
 c     CALL moyglo_pondaire(klon, zx_tmp_fi2d, airephy, 
 c    .                     ok_msk, msk, moyglo)
Index: LMDZ4/trunk/libf/phylmd/write_histhf.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_histhf.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/write_histhf.h	(revision 996)
@@ -24,7 +24,7 @@
       CALL histwrite_phy(nid_hf,"contfracATM",itau_w,zx_tmp_fi2d)
 c
-cym      CALL gr_fi_ecrit(1,klon,iim,jjmp1,pctsrf_new(:,is_ter),zx_tmp_2d)
+cym      CALL gr_fi_ecrit(1,klon,iim,jjmp1,pctsrf(:,is_ter),zx_tmp_2d)
       CALL histwrite_phy(nid_hf,"contfracOR",itau_w,
-     .                   pctsrf_new(:,is_ter))
+     .                   pctsrf(:,is_ter))
 c
 cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zt2m,zx_tmp_2d)
Index: LMDZ4/trunk/libf/phylmd/write_paramLMDZ_phy.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_paramLMDZ_phy.h	(revision 993)
+++ LMDZ4/trunk/libf/phylmd/write_paramLMDZ_phy.h	(revision 996)
@@ -7,20 +7,12 @@
 c Variables type caractere : plusieurs valeurs possibles
 c
-      IF(ocean.EQ.'force ') THEN
-       zx_tmp_2d(1:iim,1:jjmp1)=1.
-      ELSE IF(ocean.EQ.'slab  ') THEN
+      IF(type_ocean.EQ.'force ') THEN
+       zx_tmp_2d(1:iim,1:jjmp1)=1.
+      ELSE IF(type_ocean.EQ.'slab  ') THEN
        zx_tmp_2d(1:iim,1:jjmp1)=2.
-      ELSE IF(ocean.EQ.'couple') THEN
+      ELSE IF(type_ocean.EQ.'couple') THEN
        zx_tmp_2d(1:iim,1:jjmp1)=3.
       ENDIF
       CALL histwrite(nid_ctesGCM,"ocean",itau_w,
-     .               zx_tmp_2d,iim*jjmp1,ndex2d)
-c
-      IF(ok_slab_sicOBS) THEN
-       zx_tmp_2d(1:iim,1:jjmp1)=1.
-      ELSE
-       zx_tmp_2d(1:iim,1:jjmp1)=0.
-      ENDIF
-      CALL histwrite(nid_ctesGCM,"ok_slab_sicOBS",itau_w,
      .               zx_tmp_2d,iim*jjmp1,ndex2d)
 c
