Index: LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90	(revision 2073)
@@ -53,5 +53,5 @@
     ! type_ocean:      type d'ocean (force, slab, couple)
     ! version_ocean:   version d'ocean (opa8/nemo pour type_ocean=couple ou 
-    !                                   sicOBS pour type_ocean=slab)
+    !                                   sicOBS,sicINT,sicNO pour type_ocean=slab)
     ! ok_veget:   type de modele de vegetation
     ! ok_journe:  sorties journalieres
@@ -1933,5 +1933,6 @@
     IF (type_ocean=='slab' .AND. version_ocean=='xxxxxx') THEN
        version_ocean='sicOBS'
-    ELSE IF (type_ocean=='slab' .AND. version_ocean/='sicOBS') THEN
+    ELSE IF (type_ocean=='slab' .AND. version_ocean/='sicOBS' &
+             .AND. version_ocean/='sicINT' .AND. version_ocean/='sicNO') THEN
        WRITE(lunout,*)' ERROR version_ocean=',version_ocean,' not valid with slab ocean'
        CALL abort_gcm('conf_phys','version_ocean not valid',1)
Index: LMDZ5/branches/testing/libf/phylmd/dimphy.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/dimphy.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/dimphy.F90	(revision 2073)
@@ -9,6 +9,7 @@
   INTEGER,SAVE :: klevm1
   INTEGER,SAVE :: kflev
+  INTEGER,SAVE :: nslay
 
-!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon)
+!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon,nslay)
   REAL,save,allocatable,dimension(:) :: zmasq
 !$OMP THREADPRIVATE(zmasq)   
@@ -23,5 +24,5 @@
     
     klon=klon0
-    
+    nslay=1 ! Slab, provisoire (F. Codron)
     kdlon=klon
     kidia=1
Index: LMDZ5/branches/testing/libf/phylmd/limit_slab.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/limit_slab.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/limit_slab.F90	(revision 2073)
@@ -1,5 +1,5 @@
 ! $Header$
 
-SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, lmt_foce, diff_sst)
+SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, diff_sst)
 
   USE dimphy
@@ -20,10 +20,10 @@
   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
+  REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst
 
 ! Locals variables with attribute SAVE
 !****************************************************************************************
-  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, foce_save
-!$OMP THREADPRIVATE(bils_save, foce_save)
+  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, diff_sst_save
+!$OMP THREADPRIVATE(bils_save, diff_sst_save)
 
 ! Locals variables
@@ -32,6 +32,7 @@
   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
+  REAL, DIMENSION(klon_glo):: bils_glo, sst_l_glo, sst_lp1_glo, diff_sst_glo
   CHARACTER (len = 20)     :: modname = 'limit_slab'
+  LOGICAL                  :: read_bils,read_sst
 
 ! End declaration
@@ -41,12 +42,17 @@
   lmt_pas = NINT(86400./dtime)
   
+  ! F. Codron 5/14: add defaults for bils, diff_sst (0)
   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
+        read_bils=.TRUE.
+        read_sst=.TRUE.
         
         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)
+        IF (ierr /= NF90_NOERR) THEN
+            read_bils=.FALSE.
+            read_sst=.FALSE.
+        ELSE ! read file
         
         ! La tranche de donnees a lire:
@@ -57,5 +63,5 @@
 
 !****************************************************************************************
-! 2) Read bils and ocean fraction
+! 2) Read bils and SST tendency
 !
 !****************************************************************************************
@@ -63,41 +69,41 @@
 ! 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)
-!
+        IF (ierr /= NF90_NOERR) THEN
+            read_bils=.FALSE.
+        ELSE
+            ierr = NF90_GET_VAR(nid,nvarid,bils_glo,start,epais)
+            IF (ierr /= NF90_NOERR) read_bils=.FALSE.
+        END IF
 ! 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)
-
+        IF (ierr /= NF90_NOERR)  THEN
+            read_sst=.FALSE.
+        ELSE
+            ierr = NF90_GET_VAR(nid,nvarid,sst_l_glo,start,epais)
+            IF (ierr /= NF90_NOERR) read_sst=.FALSE.
 ! 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)
+            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) read_sst=.FALSE.
+        END IF
 
 !****************************************************************************************
-! 5) Close file and distribuate variables to all processus
+! 5) Close file and distribute variables to all processus
 !
 !****************************************************************************************
         ierr = NF90_CLOSE(nid)
         IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Pb when closing file', 1)
+        END IF ! Read File
+        IF (read_sst) THEN
+! 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)
+            DO i=1, klon_glo
+               diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
+            END DO
+        END IF !read_sst
      ENDIF ! is_mpi_root
 
@@ -105,18 +111,24 @@
        
      IF (.NOT. ALLOCATED(bils_save)) THEN
-        ALLOCATE(bils_save(klon), foce_save(klon), stat=ierr)
+        ALLOCATE(bils_save(klon), diff_sst_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)
+! Giveddefault values if needed
+     IF (read_bils) THEN
+         CALL Scatter(bils_glo, bils_save)
+     ELSE
+         bils_save(:)=0.
+     END IF
+     IF (read_sst) THEN
+         CALL Scatter(diff_sst_glo, diff_sst_save)
+     ELSE
+         diff_sst_save(:)=0.
+     END IF
      
-  ELSE ! not time to read
-     diff_sst(:) = 0.
   ENDIF ! time to read
 
   lmt_bils(:) = bils_save(:)
-  lmt_foce(:) = foce_save(:)
+  diff_sst(:) = diff_sst_save(:)
   
 END SUBROUTINE limit_slab
Index: LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90	(revision 2073)
@@ -5,7 +5,26 @@
 ! "ocean=slab".
 !
+
+  USE dimphy
+  USE indice_sol_mod
+
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: ocean_slab_frac, ocean_slab_noice
+  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice!, ocean_slab_ice
+
+  INTEGER, PRIVATE, SAVE                           :: cpl_pas
+  !$OMP THREADPRIVATE(cpl_pas)
+  REAL, PRIVATE, SAVE                              :: cyang
+  !$OMP THREADPRIVATE(cyang)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
+  !$OMP THREADPRIVATE(slabh)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
+  !$OMP THREADPRIVATE(tslab)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE  :: pctsrf
+  !$OMP THREADPRIVATE(pctsrf)
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
+  !$OMP THREADPRIVATE(slab_bils)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
+  !$OMP THREADPRIVATE(bils_cum)
 
 CONTAINS
@@ -13,10 +32,74 @@
 !****************************************************************************************
 !
-  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
-
-    USE dimphy
+  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
+  !, seaice_rst etc
+
+    use IOIPSL
+
+    INCLUDE "iniprint.h"
+    ! For ok_xxx vars (Ekman...)
+    INCLUDE "clesphys.h"
+
+    ! Input variables
+!****************************************************************************************
+    REAL, INTENT(IN)                         :: dtime
+! Variables read from restart file
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
+
+! Local variables
+!****************************************************************************************
+    INTEGER                :: error
+    CHARACTER (len = 80)   :: abort_message
+    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
+
+!****************************************************************************************
+! Allocate surface fraction read from restart file
+!****************************************************************************************
+    ALLOCATE(pctsrf(klon,nbsrf), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_pctsrf_slab'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    pctsrf(:,:) = pctsrf_rst(:,:)
+
+!****************************************************************************************
+! Allocate local variables
+!****************************************************************************************
+    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(bils_cum(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slab_bils_cum'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    bils_cum(:) = 0.0   
+
+! Layer thickness
+    ALLOCATE(slabh(nslay), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slabh'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    slabh(1)=50.
+! cyang = 1/heat capacity of top layer (rho.c.H)
+    cyang=1/(slabh(1)*4.228e+06)
+
+! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
+! pour un calcul à chaque pas de temps, cpl_pas=1
+    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
+    CALL getin('cpl_pas',cpl_pas)
+    print *,'cpl_pas',cpl_pas
+ END SUBROUTINE ocean_slab_init
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
+
     USE limit_read_mod
     USE surface_data
-    USE indice_sol_mod
 
 !    INCLUDE "clesphys.h"
@@ -27,20 +110,15 @@
     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
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
     LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
 
 ! Local variables
 !****************************************************************************************
-    CHARACTER (len = 80)   :: abort_message
-    CHARACTER (len = 20)   :: modname = 'ocean_slab_frac'
-
-
-    IF (version_ocean == 'sicOBS') THEN   
-       CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
+
+    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
+       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, 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).  
+       pctsrf_chg(:,:)=pctsrf(:,:)
+       is_modified=.TRUE.
     END IF
 
@@ -57,9 +135,8 @@
        radsol, snow, agesno, &
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, lmt_bils)
+       tsurf_new, dflux_s, dflux_l, qflux)
     
-    USE dimphy
     USE calcul_fluxs_mod
-    USE indice_sol_mod
+    USE surface_data
 
     INCLUDE "iniprint.h"
@@ -95,15 +172,13 @@
     REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)   :: lmt_bils
+    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
 
 ! Local variables
 !****************************************************************************************
-    INTEGER               :: i
+    INTEGER               :: i,ki
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd
-    REAL, DIMENSION(klon) :: lmt_bils_oce, lmt_foce, diff_sst
+    REAL, DIMENSION(klon) :: diff_sst, lmt_bils
     REAL, DIMENSION(klon) :: u0, v0
     REAL, DIMENSION(klon) :: u1_lay, v1_lay
-    REAL                  :: calc_bils_oce, deltat
-    REAL, PARAMETER       :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
 
 !****************************************************************************************
@@ -136,28 +211,124 @@
          flux_u1, flux_v1)  
 
-!****************************************************************************************
-! 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
+! Accumulate total fluxes locally
+    slab_bils(:)=0.
+    DO i=1,knon
+        ki=knindex(i)
+        slab_bils(ki)=(fluxlat(i)+fluxsens(i)+radsol(i))*pctsrf(ki,is_oce)/(1.-zmasq(ki))
+        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
+! Also taux, tauy, saved vars...
+    END DO
+
+!****************************************************************************************
+! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
+!
+!****************************************************************************************
+    lmt_bils(:)=0.
+    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst)  ! global pour un processus
+    ! lmt_bils and diff_sst saved by limit_slab
+    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.
+    ! qflux = total QFlux correction (in W/m2)
 
 !****************************************************************************************
 ! 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
-
+!***********************************************o*****************************************
+    tsurf_new=tsurf_in
+    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
+        ! Compute transport
+        ! Add read QFlux and SST tendency
+        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 
+        ! Add cumulated surface fluxes
+        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
+        ! Update surface temperature
+        SELECT CASE(version_ocean)
+        CASE('sicNO')
+            DO i=1,knon
+                ki=knindex(i)
+                tsurf_new(i)=tslab(ki,1)
+            END DO
+        CASE('sicOBS') ! check for sea ice or tsurf below freezing
+            DO i=1,knon
+                ki=knindex(i)
+                IF ((tslab(ki,1).LT.t_freeze).OR.(pctsrf(ki,is_sic).GT.epsfra)) THEN
+                    tsurf_new(i)=t_freeze
+                    tslab(ki,1)=t_freeze
+                ELSE
+                    tsurf_new(i)=tslab(ki,1)
+                END IF
+            END DO
+        CASE('sicINT')
+            DO i=1,knon
+                ki=knindex(i)
+                IF (pctsrf(ki,is_sic).LT.epsfra) THEN ! Free of ice
+                    IF (tslab(ki,1).GT.t_freeze) THEN 
+                        tsurf_new(i)=tslab(ki,1)
+                    ELSE
+                        tsurf_new(i)=t_freeze
+                        ! Call new ice routine
+                        tslab(ki,1)=t_freeze
+                    END IF
+                ELSE ! ice present, tslab update completed in slab_ice
+                    tsurf_new(i)=t_freeze
+                END IF !ice free
+            END DO
+        END SELECT
+        bils_cum(:)=0.0! clear cumulated fluxes
+    END IF ! coupling time
   END SUBROUTINE ocean_slab_noice
 !
 !****************************************************************************************
 !
+!  SUBROUTINE ocean_slab_ice(   &
+!       itime, dtime, jour, knon, knindex, &
+!       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+!       AcoefH, AcoefQ, BcoefH, BcoefQ, &
+!       AcoefU, AcoefV, BcoefU, BcoefV, &
+!       ps, u1, v1, &
+!       radsol, snow, qsurf, qsol, agesno, tsoil, &
+!       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
+!       tsurf_new, dflux_s, dflux_l)
+!
+!****************************************************************************************
+! 1) Flux calculation
+!****************************************************************************************
+! set beta, cal etc. depends snow / ice surf ?
+! calcul_fluxs (sens, lat etc)
+! calcul_flux_wind
+
+!****************************************************************************************
+! 2) Update surface
+!****************************************************************************************
+! neige, fonte
+! flux glace-ocean
+! update temperature
+! neige precip, evap
+! Melt snow & ice from above
+! New albedo
+
+!****************************************************************************************
+! 3) Recalculate new ocean temperature
+!    Melt / freeze from below
+!***********************************************o*****************************************
+
+
+!  END SUBROUTINE ocean_slab_ice
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_final
+  !, seaice_rst etc
+
+    ! For ok_xxx vars (Ekman...)
+    INCLUDE "clesphys.h"
+
+!****************************************************************************************
+! Deallocate module variables
+!
+!****************************************************************************************
+    IF (ALLOCATED(pctsrf)) DEALLOCATE(pctsrf)
+    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
+
+  END SUBROUTINE ocean_slab_final
+
 END MODULE ocean_slab_mod
Index: LMDZ5/branches/testing/libf/phylmd/phyetat0.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phyetat0.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/phyetat0.F90	(revision 2073)
@@ -3,5 +3,5 @@
 SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
 
-  USE dimphy, only: klon, zmasq, klev
+  USE dimphy, only: klon, zmasq, klev, nslay
   USE iophy, ONLY : init_iophy_new
   USE ocean_cpl_mod,    ONLY : ocean_cpl_init
@@ -15,5 +15,5 @@
        solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
        wake_deltat, wake_fip, wake_pe, wake_s, zgam, zmax0, zmea, zpic, zsig, &
-       zstd, zthe, zval
+       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
   USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
   USE infotrac, only: nbtr, type_trac, tname, niadv
@@ -21,4 +21,5 @@
   USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
   USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
+  USE ocean_slab_mod, ONLY: tslab, ocean_slab_init
 
   IMPLICIT none
@@ -40,5 +41,4 @@
 
   REAL tsoil(klon, nsoilmx, nbsrf)
-  REAL tslab(klon), seaice(klon)
   REAL qsurf(klon, nbsrf)
   REAL qsol(klon)
@@ -997,4 +997,40 @@
   PRINT*, '(ecart-type) detr_therm:', xmin, xmax
 
+  CALL get_field("ALE_BL", ale_bl, found)
+  IF (.NOT. found) THEN
+     PRINT*, "phyetat0: Le champ <ale_bl> est absent"
+     PRINT*, "Depart legerement fausse. Mais je continue"
+     ale_bl=0.
+  ENDIF
+  xmin = 1.0E+20
+  xmax = -1.0E+20
+  xmin = MINval(ale_bl)
+  xmax = MAXval(ale_bl)
+  PRINT*, '(ecart-type) ale_bl:', xmin, xmax
+
+  CALL get_field("ALE_BL_TRIG", ale_bl_trig, found)
+  IF (.NOT. found) THEN
+     PRINT*, "phyetat0: Le champ <ale_bl_trig> est absent"
+     PRINT*, "Depart legerement fausse. Mais je continue"
+     ale_bl_trig=0.
+  ENDIF
+  xmin = 1.0E+20
+  xmax = -1.0E+20
+  xmin = MINval(ale_bl_trig)
+  xmax = MAXval(ale_bl_trig)
+  PRINT*, '(ecart-type) ale_bl_trig:', xmin, xmax
+
+  CALL get_field("ALP_BL", alp_bl, found)
+  IF (.NOT. found) THEN
+     PRINT*, "phyetat0: Le champ <alp_bl> est absent"
+     PRINT*, "Depart legerement fausse. Mais je continue"
+     alp_bl=0.
+  ENDIF
+  xmin = 1.0E+20
+  xmax = -1.0E+20
+  xmin = MINval(alp_bl)
+  xmax = MAXval(alp_bl)
+  PRINT*, '(ecart-type) alp_bl:', xmin, xmax
+
   ! Read and send field trs to traclmdz
 
@@ -1045,4 +1081,21 @@
   end if
 
+  ! Initialize Slab variables
+  IF ( type_ocean == 'slab' ) THEN
+      ALLOCATE(tslab(klon,nslay), stat=ierr)
+      IF (ierr /= 0) CALL abort_gcm &
+         ('phyetat0', 'pb allocation tslab', 1)
+      CALL get_field("tslab", tslab, found)
+      IF (.NOT. found) THEN 
+          PRINT*, "phyetat0: Le champ <tslab> est absent"
+          PRINT*, "Initialisation a tsol_oce"
+          DO i=1,nslay
+              tslab(:,i)=ftsol(:,is_oce)
+          END DO
+      END IF 
+      print*, "calling slab_init"
+      CALL ocean_slab_init(dtime, pctsrf)
+  END IF ! Slab        
+
   ! on ferme le fichier
   CALL close_startphy
Index: LMDZ5/branches/testing/libf/phylmd/phyredem.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phyredem.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/phyredem.F90	(revision 2073)
@@ -15,4 +15,6 @@
   USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
   USE indice_sol_mod
+  USE surface_data
+  USE ocean_slab_mod, ONLY : tslab
 
   IMPLICIT none
@@ -33,5 +35,4 @@
 
   REAL tsoil(klon, nsoilmx, nbsrf)
-  REAL tslab(klon), seaice(klon)
   REAL qsurf(klon, nbsrf)
   REAL qsol(klon)
@@ -322,4 +323,10 @@
 
   CALL put_field("DETR_THERM", "DETR_THERM", detr_therm)
+
+  CALL put_field("ALE_BL", "ALE_BL", Ale_bl)
+
+  CALL put_field("ALE_BL_TRIG", "ALE_BL_TRIG", Ale_bl_trig)
+
+  CALL put_field("ALP_BL", "ALP_BL", Alp_bl)
 
   ! trs from traclmdz_mod
@@ -340,4 +347,9 @@
   END IF
 
+  ! Restart variables for Slab ocean
+  IF (type_ocean == 'slab') THEN
+      CALL put_field("tslab", "Slab ocean temperature", tslab)
+  END IF
+
   if (ok_gwd_rando) then
      call put_field("du_gwd_rando", &
Index: LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2073)
@@ -456,5 +456,9 @@
     !(/ ('', i=1, 9) /))
   TYPE(ctrl_out), SAVE :: o_slab_bils = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
-    'slab_bils_oce', 'Bilan au sol sur ocean slab', 'W/m2', (/ ('', i=1, 9) /))
+    'slab_bils', 'flux atmos - slab ponderes foce', 'W/m2', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_qflux = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'slab_qflux', 'Correction flux slab', 'W/m2', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_tslab = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'tslab', 'Temperature ocean slab', 'K', (/ ('', i=1, 9) /))
   TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), &
     'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /))
Index: LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90	(revision 2073)
@@ -24,5 +24,5 @@
     ! defined and initialised in phys_output_mod.F90
 
-    USE dimphy, only: klon, klev, klevp1
+    USE dimphy, only: klon, klev, klevp1, nslay
     USE control_mod, only: day_step, iphysiq
     USE phys_output_ctrlout_mod, only: o_phis, o_aire, is_ter, is_lic, is_oce, &
@@ -77,5 +77,6 @@
          o_ale_bl_trig, o_alp_bl_det, &
          o_alp_bl_fluct_m, o_alp_bl_fluct_tke, &
-         o_alp_bl_conv, o_alp_bl_stat, o_slab_bils, &
+         o_alp_bl_conv, o_alp_bl_stat, &
+         o_slab_qflux, o_tslab, o_slab_bils, &
          o_weakinv, o_dthmin, o_cldtau, &
          o_cldemi, o_pr_con_l, o_pr_con_i, &
@@ -219,4 +220,5 @@
          bils_ec,bils_ech, bils_tke, bils_kinetic, bils_latent, bils_enthalp, &
          itau_con, nfiles, clef_files, nid_files, zvstr_gwd_rando
+    USE ocean_slab_mod, only: tslab, slab_bils
     USE indice_sol_mod, only: nbsrf
     USE infotrac, only: nqtot
@@ -726,6 +728,14 @@
        ENDIF  !(iflag_clos_bl>=1)
 !!! fin nrlmd le 10/04/2012
+       ! Output of slab ocean variables
        IF (type_ocean=='slab ') THEN
-          CALL histwrite_phy(o_slab_bils, slab_wfbils)
+          CALL histwrite_phy(o_slab_qflux, slab_wfbils)
+          CALL histwrite_phy(o_slab_bils, slab_bils)
+          IF (nslay.EQ.1) THEN
+              zx_tmp_fi2d(:)=tslab(:,1)
+              CALL histwrite_phy(o_tslab, zx_tmp_fi2d)
+          ELSE
+              CALL histwrite_phy(o_tslab, tslab)
+          END IF
        ENDIF !type_ocean == force/slab
        CALL histwrite_phy(o_weakinv, weak_inversion)
Index: LMDZ5/branches/testing/libf/phylmd/physiq.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/physiq.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/physiq.F90	(revision 2073)
@@ -442,5 +442,5 @@
   REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
   ! RomP <<<
-  
+
   REAL          :: calday
 
@@ -489,6 +489,6 @@
   EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
   !AA 
-! JBM (3/14) fisrtilp_tr not loaded
-! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
+  ! JBM (3/14) fisrtilp_tr not loaded
+  ! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
   !                          ! stockage des coefficients necessaires au
   !                          ! lessivage OFF-LINE et ON-LINE
@@ -542,5 +542,5 @@
 
   !
-!  REAL zxsnow(klon)
+  !  REAL zxsnow(klon)
   REAL zxsnow_dummy(klon)
   !
@@ -636,5 +636,5 @@
   !======================================================================
   !
-  
+
   !
   integer itau_w   ! pas de temps ecriture = itap + itau_phy
@@ -689,7 +689,7 @@
   !
   !IM 280405 BEG
-!  INTEGER nid_bilKPins, nid_bilKPave
-!  SAVE nid_bilKPins, nid_bilKPave
-!  !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
+  !  INTEGER nid_bilKPins, nid_bilKPave
+  !  SAVE nid_bilKPins, nid_bilKPave
+  !  !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
   !
   REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
@@ -701,5 +701,5 @@
   REAL zsto
   REAL zstophy, zout
- 
+
   real zjulian
   save zjulian
@@ -1121,8 +1121,8 @@
         endif
 
-        do i = 1,klon
-           Ale_bl(i)=0.
-           Alp_bl(i)=0.
-        enddo
+!        do i = 1,klon
+!           Ale_bl(i)=0.
+!           Alp_bl(i)=0.
+!        enddo
 
         !================================================================================
@@ -1588,8 +1588,8 @@
 
      IF (solarlong0<-999.) then
-     ! Generic case with evolvoing season
+        ! Generic case with evolvoing season
         zzz=real(days_elapsed+1)
      ELSE IF (abs(solarlong0-1000.)<1.e-4) then
-     ! Particular case with annual mean insolation
+        ! Particular case with annual mean insolation
         zzz=real(90) ! could be revisited
         IF (read_climoz/=-1) THEN
@@ -1598,5 +1598,5 @@
         ENDIF
      ELSE
-     ! Case where the season is imposed with solarlong0
+        ! Case where the season is imposed with solarlong0
         zzz=real(90) ! could be revisited
      ENDIF
@@ -1764,8 +1764,8 @@
      IF (klon_glo==1) THEN
         CALL add_pbl_tend &
-          (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
+             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
      ELSE
         CALL add_phys_tend &
-          (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
+             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
      ENDIF
      !--------------------------------------------------------------------
@@ -1987,4 +1987,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+
      endif
      do i=1,klon
@@ -2369,44 +2370,44 @@
                    s_trig,s2,n2
            ENDIF
- 
-!Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
+
+           !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
            IF (iflag_trig_bl.eq.1) then
 
-           !----Tirage al\'eatoire et calcul de ale_bl_trig
-           do i=1,klon
-              if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
-                 proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
-                      (n2(i)*dtime/tau_trig(i))
-                 !        print *, 'proba_notrig(i) ',proba_notrig(i)
-                 if (random_notrig(i) .ge. proba_notrig(i)) then 
-                    ale_bl_trig(i)=ale_bl_stat(i)
+              !----Tirage al\'eatoire et calcul de ale_bl_trig
+              do i=1,klon
+                 if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
+                    proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
+                         (n2(i)*dtime/tau_trig(i))
+                    !        print *, 'proba_notrig(i) ',proba_notrig(i)
+                    if (random_notrig(i) .ge. proba_notrig(i)) then 
+                       ale_bl_trig(i)=ale_bl_stat(i)
+                    else
+                       ale_bl_trig(i)=0.
+                    endif
                  else
+                    proba_notrig(i)=1.
+                    random_notrig(i)=0.
                     ale_bl_trig(i)=0.
                  endif
-              else
-                 proba_notrig(i)=1.
-                 random_notrig(i)=0.
-                 ale_bl_trig(i)=0.
-              endif
-           enddo
+              enddo
 
            ELSE IF (iflag_trig_bl.eq.2) then
 
-           do i=1,klon
-              if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
-                 proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
-                      (n2(i)*dtime/tau_trig(i))
-                 !        print *, 'proba_notrig(i) ',proba_notrig(i)
-                 if (random_notrig(i) .ge. proba_notrig(i)) then 
-                    ale_bl_trig(i)=Ale_bl(i)
+              do i=1,klon
+                 if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
+                    proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
+                         (n2(i)*dtime/tau_trig(i))
+                    !        print *, 'proba_notrig(i) ',proba_notrig(i)
+                    if (random_notrig(i) .ge. proba_notrig(i)) then 
+                       ale_bl_trig(i)=Ale_bl(i)
+                    else
+                       ale_bl_trig(i)=0.
+                    endif
                  else
+                    proba_notrig(i)=1.
+                    random_notrig(i)=0.
                     ale_bl_trig(i)=0.
                  endif
-              else
-                 proba_notrig(i)=1.
-                 random_notrig(i)=0.
-                 ale_bl_trig(i)=0.
-              endif
-           enddo
+              enddo
 
            ENDIF
@@ -2424,13 +2425,13 @@
 
            do i=1,klon
-!CR: alp probabiliste
-               if (ale_bl_trig(i).gt.0.) then
-                  alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
-               endif
-           enddo        
-  
+              !CR: alp probabiliste
+              if (ale_bl_trig(i).gt.0.) then
+                 alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
+              endif
+           enddo
+
         else if (iflag_clos_bl.eq.2) then
 
-!CR: alp calculee dans thermcell_main
+           !CR: alp calculee dans thermcell_main
            do i=1,klon
               alp_bl(i)=alp_bl_stat(i)
@@ -2469,7 +2470,7 @@
 
         do i=1,klon
-!           zmax_th(i)=pphi(i,lmax_th(i))/rg
-!CR:04/05/12:correction calcul zmax
-         zmax_th(i)=zmax0(i) 
+           !           zmax_th(i)=pphi(i,lmax_th(i))/rg
+           !CR:04/05/12:correction calcul zmax
+           zmax_th(i)=zmax0(i) 
         enddo
 
@@ -2500,9 +2501,9 @@
 
         if (iflag_thermals==0) then
-        ! Calling adjustment alone (but not the thermal plume model)
+           ! Calling adjustment alone (but not the thermal plume model)
            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
                 , d_t_ajsb, d_q_ajsb)
         else if (iflag_thermals>0) then
-        ! Calling adjustment above the top of thermal plumes
+           ! Calling adjustment above the top of thermal plumes
            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
                 , d_t_ajsb, d_q_ajsb)
@@ -2678,27 +2679,27 @@
         IF (.NOT. aerosol_couple) THEN
            IF (iflag_rrtm .EQ. 0) THEN !--old radiation
-!
-             CALL readaerosol_optic( &
-             debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
-             pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
-             mass_solu_aero, mass_solu_aero_pi,  &
-             tau_aero, piz_aero, cg_aero,  &
-             tausum_aero, tau3d_aero)
-!
+              !
+              CALL readaerosol_optic( &
+                   debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
+                   pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
+                   mass_solu_aero, mass_solu_aero_pi,  &
+                   tau_aero, piz_aero, cg_aero,  &
+                   tausum_aero, tau3d_aero)
+              !
            ELSE                       ! RRTM radiation
-!
+              !
 #ifdef CPP_RRTM
-             CALL readaerosol_optic_rrtm( &
-             debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
-             pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
-             mass_solu_aero, mass_solu_aero_pi,  &
-             tau_aero_rrtm, piz_aero_rrtm, cg_aero_rrtm,  &
-             tausum_aero, tau3d_aero)
+              CALL readaerosol_optic_rrtm( &
+                   debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
+                   pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
+                   mass_solu_aero, mass_solu_aero_pi,  &
+                   tau_aero_rrtm, piz_aero_rrtm, cg_aero_rrtm,  &
+                   tausum_aero, tau3d_aero)
 #else
 
-            abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
-            call abort_gcm(modname,abort_message,1)
+              abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
+              call abort_gcm(modname,abort_message,1)
 #endif
-!
+              !
            ENDIF
         ENDIF
@@ -2706,11 +2707,11 @@
         tausum_aero(:,:,:) = 0.
         IF (iflag_rrtm .EQ. 0) THEN !--old radiation
-          tau_aero(:,:,:,:) = 0.
-          piz_aero(:,:,:,:) = 0.
-          cg_aero(:,:,:,:)  = 0.
+           tau_aero(:,:,:,:) = 0.
+           piz_aero(:,:,:,:) = 0.
+           cg_aero(:,:,:,:)  = 0.
         ELSE
-          tau_aero_rrtm(:,:,:,:)=0.0
-          piz_aero_rrtm(:,:,:,:)=0.0
-          cg_aero_rrtm(:,:,:,:)=0.0
+           tau_aero_rrtm(:,:,:,:)=0.0
+           piz_aero_rrtm(:,:,:,:)=0.0
+           cg_aero_rrtm(:,:,:,:)=0.0
         ENDIF
      ENDIF
@@ -2721,12 +2722,12 @@
         PRINT *,'appel a readaerosolstrat', mth_cur
         IF (iflag_rrtm.EQ.0) THEN
-         CALL readaerosolstrato(debut)
+           CALL readaerosolstrato(debut)
         ELSE
 #ifdef CPP_RRTM
-         CALL readaerosolstrato_rrtm(debut)
+           CALL readaerosolstrato_rrtm(debut)
 #else
 
-         abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
-         call abort_gcm(modname,abort_message,1)
+           abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
+           call abort_gcm(modname,abort_message,1)
 #endif
         ENDIF
@@ -2951,12 +2952,12 @@
      IF (iflag_rrtm.NE.0) THEN 
 #ifdef CPP_RRTM
-       IF (ok_cdnc.AND.NRADLP.NE.3) THEN
-         abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc' 
-         call abort_gcm(modname,abort_message,1)
-       endif
+        IF (ok_cdnc.AND.NRADLP.NE.3) THEN
+           abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc' 
+           call abort_gcm(modname,abort_message,1)
+        endif
 #else
 
-       abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
-       call abort_gcm(modname,abort_message,1)
+        abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
+        call abort_gcm(modname,abort_message,1)
 #endif
      ENDIF
Index: LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 2073)
@@ -12,4 +12,5 @@
   USE aero_mod
   USE phys_local_var_mod, only: od550aer,od865aer,ec550aer,od550lt1aer
+  USE YOMCST           , only : RD , RG
 
   !
@@ -49,5 +50,4 @@
   !
   IMPLICIT NONE
-  INCLUDE "YOMCST.h"
   !
   ! Input arguments:
Index: LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 2073)
@@ -11,4 +11,5 @@
   USE aero_mod
   USE phys_local_var_mod, only: absvisaer
+  USE YOMCST            , only: RD , RG
 
   !    Yves Balkanski le 12 avril 2006
@@ -21,5 +22,4 @@
   IMPLICIT NONE
 
-  INCLUDE "YOMCST.h"
   INCLUDE "iniprint.h"
   INCLUDE "clesphys.h"
Index: LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato_rrtm.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato_rrtm.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato_rrtm.F90	(revision 2073)
@@ -19,5 +19,4 @@
     implicit none
 
-    include "YOMCST.h"
     include "dimensions.h"
 
Index: LMDZ5/branches/testing/libf/phylmd/rrtm/yoesw.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/rrtm/yoesw.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/rrtm/yoesw.F90	(revision 2073)
@@ -15,4 +15,5 @@
 REAL(KIND=JPRB) :: RRAY(6,6)
 REAL(KIND=JPRB), ALLOCATABLE :: RSUN(:)
+!$OMP THREADPRIVATE(RSUN)
 REAL(KIND=JPRB) :: RPDH1
 REAL(KIND=JPRB) :: RPDU1
Index: LMDZ5/branches/testing/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/surf_ocean_mod.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/surf_ocean_mod.F90	(revision 2073)
@@ -138,4 +138,15 @@
 
 !****************************************************************************************
+! fcodron: compute lmt_bils  forced case (same as wfbils_oce / 1.-contfracatm)
+!****************************************************************************************
+    IF (type_ocean.NE.'slab') THEN
+        lmt_bils(:)=0.
+        DO i=1,knon
+           lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
+           *pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i))) 
+        END DO
+    END IF
+
+!****************************************************************************************
 ! Calculate albedo
 !
Index: LMDZ5/branches/testing/libf/phylmd/surf_seaice_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/surf_seaice_mod.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/surf_seaice_mod.F90	(revision 2073)
@@ -108,5 +108,16 @@
             tsurf_new, dflux_s, dflux_l)
        
-    ELSE IF (type_ocean == 'force' .OR. (type_ocean == 'slab' .AND. version_ocean=='sicOBS')) THEN
+!    ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
+!       CALL ocean_slab_ice( & 
+!          itime, dtime, jour, knon, knindex, &
+!          debut, &
+!          tsurf, p1lay, cdragh, precip_rain, precip_snow, temp_air, spechum,&
+!          AcoefH, AcoefQ, BcoefH, BcoefQ, &
+!          ps, u1, v1, &
+!          radsol, snow, qsurf, qsol, agesno, tsoil, &
+!          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+!          tsurf_new, dflux_s, dflux_l)
+!
+      ELSE ! type_ocean=force or slab +sicOBS or sicNO
        CALL ocean_forced_ice( &
             itime, dtime, jour, knon, knindex, &
@@ -118,15 +129,4 @@
             qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
             tsurf_new, dflux_s, dflux_l)
-
-    ELSE IF (type_ocean == 'slab') THEN
-!!$       CALL ocean_slab_ice( & 
-!!$          itime, dtime, jour, knon, knindex, &
-!!$          debut, &
-!!$          tsurf, p1lay, cdragh, precip_rain, precip_snow, temp_air, spechum,&
-!!$          AcoefH, AcoefQ, BcoefH, BcoefQ, &
-!!$          ps, u1, v1, pctsrf, &
-!!$          radsol, snow, qsurf, qsol, agesno, tsoil, &
-!!$          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-!!$          tsurf_new, dflux_s, dflux_l)
 
     END IF
Index: LMDZ5/branches/testing/libf/phylmd/surface_data.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/surface_data.F90	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/surface_data.F90	(revision 2073)
@@ -7,4 +7,6 @@
   REAL, PARAMETER        :: tau_gl=86400.*5.
   REAL, PARAMETER        :: calsno=1./(2.3867e+06*.15)
+  REAL, PARAMETER        :: t_freeze=271.35
+  REAL, PARAMETER        :: t_melt=273.15
   
   LOGICAL, SAVE          :: ok_veget      ! true for use of vegetation model ORCHIDEE
@@ -20,5 +22,5 @@
 
   ! if type_ocean=couple : version_ocean=opa8 ou nemo
-  ! if type_ocean=slab   : version_ocean=sicOBS
+  ! if type_ocean=slab   : version_ocean=sicOBS or sicINT or sicNO
   CHARACTER(len=6), SAVE :: version_ocean 
   !$OMP THREADPRIVATE(version_ocean)
