Index: trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90	(revision 3342)
+++ trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90	(revision 3352)
@@ -106,4 +106,5 @@
     REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice [in W/(m.K) or (W.kg)/(K.m4)]
     REAL, PRIVATE, SAVE :: sno_cond ! conductivity of snow [in W/(m.K) or (W.kg)/(K.m4)]
+!$OMP THREADPRIVATE(sno_cond)
     REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice [in J/(kg.K)]
     REAL, PARAMETER :: sea_cap=3994.   ! specific heat capacity, seawater [in J/(kg.K)]
@@ -116,4 +117,5 @@
     REAL, PARAMETER :: ice_frac_min=0.001
     REAL, PRIVATE, SAVE :: ice_frac_max ! Max ice fraction (leads)
+!$OMP THREADPRIVATE(ice_frac_max)
     REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness [in kg/m2]
     REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness [in kg/m2]
@@ -121,4 +123,5 @@
     ! ice_thin is also height of new ice
     REAL, PRIVATE, SAVE :: h_ice_thick ! thin ice thickness
+!$OMP THREADPRIVATE(h_ice_thick)
     ! above ice_thick, priority is melt height / grow lateral
     REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice [in kg/m2]
@@ -1049,6 +1052,7 @@
 ! ------ End Albedo ----------
 
-    !tests remaining ice fraction
-    WHERE ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min))
+    !tests remaining ice fraction (on ocean grid points)
+    WHERE ((zmasq==0).and. &
+           ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min)))
         tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
         tice=t_melt
@@ -1085,6 +1089,6 @@
 !****************************************************************************************
 !
-  SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_g_loc, &
-       dt_hdiff_loc,dt_ekman_loc)
+  SUBROUTINE ocean_slab_get_vars(ngrid, tslab_loc, tice_loc, seaice_loc, &
+                                 flux_g_loc, dt_hdiff_loc,dt_ekman_loc)
 
 ! "Get some variables from module ocean_slab_mod"
@@ -1092,4 +1096,5 @@
     INTEGER, INTENT(IN)                     :: ngrid
     REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: tslab_loc
+    REAL, INTENT(OUT) :: tice_loc(ngrid)
     REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc
     REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc
@@ -1100,10 +1105,10 @@
 
 ! Set the output variables
-    tslab_loc(:,:) = 0.
+    tslab_loc(:,:) = tslab(:,:)
+    tice_loc(:)=tice(:)
+    seaice_loc(:) = seaice(:)
+    flux_g_loc(:) = slab_bilg(:)
     dt_hdiff_loc(:,:)=0.
     dt_ekman_loc(:,:)=0.
-    tslab_loc(:,:) = tslab(:,:)
-    seaice_loc(:) = seaice(:)
-    flux_g_loc(:) = slab_bilg(:)
 
 !!      dt_hdiff_loc(:,i) = dt_hdiff(:,i)*slabh(i)*1000.*4228. !Convert en W/m2
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3342)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3352)
@@ -666,5 +666,8 @@
 
 !!           CALL init_masquv(ngrid,zmasq)
-
+           ! ensure that "slab" variables match "physiq" ones
+           call ocean_slab_get_vars(ngrid, tslab, tsea_ice, sea_ice, flux_g, &
+                      dt_hdiff, dt_ekman)
+           
          endif ! end of 'ok_slab_ocean'.
 
@@ -678,15 +681,23 @@
 
             if (ok_slab_ocean) then
-               do ig=1,ngrid
-                  if (nint(rnat(ig)).eq.2) then
-                     capcal(ig)=capcalocean
-                     if (pctsrf_sic(ig).gt.0.5) then
-                        capcal(ig)=capcalseaice
-                        if (qsurf(ig,igcm_h2o_ice).gt.0.) then
-                           capcal(ig)=capcalsno
-                        endif
-                     endif
-                  endif
-               enddo
+              ! mimic what is done during a regular time step
+              do ig=1,ngrid
+                fluxgrdocean(ig)=fluxgrd(ig)
+                if (nint(rnat(ig)).eq.0) then
+                  capcal(ig)=capcalocean
+                  fluxgrd(ig)=0.
+                  ! Dividing by cell area to have flux in W/m2
+                  fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))/cell_area(ig)
+                  do iq=1,nsoilmx
+                    tsoil(ig,iq)=tsurf(ig)
+                  enddo
+                  if (pctsrf_sic(ig).gt.0.5) then
+                    capcal(ig)=capcalseaice
+                    if (qsurf(ig,igcm_h2o_ice).gt.0.) then
+                      capcal(ig)=capcalsno
+                    endif
+                  endif ! of if (pctsrf_sic(ig).gt.0.5)
+                endif ! of if (nint(rnat(ig)).eq.0)
+              enddo ! of do ig=1,ngrid
             endif ! end of 'ok_slab_ocean'.
 
@@ -1979,5 +1990,5 @@
             call ocean_slab_frac(pctsrf_sic, zmasq)
 
-            call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, &
+            call ocean_slab_get_vars(ngrid, tslab, tsea_ice, sea_ice, flux_g, &
                       dt_hdiff, dt_ekman)
 		      
@@ -2374,4 +2385,8 @@
          if (ngrid.ne.1) then
             write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+            
+            ! fetch "ocean variables" to ensure they are stored
+            call ocean_slab_get_vars(ngrid, tslab, tsea_ice, sea_ice, flux_g, &
+                                     dt_hdiff, dt_ekman)
 
             call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
