Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5115)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5116)
@@ -184,9 +184,17 @@
             ENDIF
 
-            !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
+            ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
             zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
 
+            ! --- Energy from the turbulent oceanic heat flux (W/m2) --- !
+            zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 
+            fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
+            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
+            ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach 
+            !                        the freezing point, so that we do not have SST < T_freeze
+            !                        This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
+
             !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice
-            qlead(ji,jj) = MIN( 0._wp , zqld - zqfr ) 
+            qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
 
             ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 
@@ -198,13 +206,4 @@
             ENDIF
             !
-            !-- Energy from the turbulent oceanic heat flux --- !
-            !clem zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )
-            zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 
-            fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 
-            ! upper bound for fhtur: we do not want SST to drop below Tfreeze. 
-            ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr)   
-            ! This is not a clean budget, so that should be corrected at some point
-            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
-
             ! -----------------------------------------
             ! Net heat flux on top of ice-ocean [W.m-2]
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5115)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5116)
@@ -91,5 +91,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
       REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
-      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt     ! Heat budget to determine melting or freezing(W.m-2)
+      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
       INTEGER , POINTER, DIMENSION(:) ::   icount      ! number of layers vanished by melting 
 
@@ -500,5 +500,5 @@
       DO jk = nlay_i, 1, -1
          DO ji = kideb, kiut
-            IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared from surface melting 
+            IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared by surface melting 
 
                ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5115)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5116)
@@ -480,5 +480,5 @@
                rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
                !
-               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1.0 - rswitch )
+               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
                !
                zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
@@ -500,5 +500,5 @@
          DO jk = 1, nlay_i
             zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
-            zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
+            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
             DO ji = kideb, kiut
                s_i_1d(ji,jk) = zsal
@@ -516,12 +516,15 @@
       !!                   ***  ROUTINE lim_var_zapsmall ***
       !!
-      !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes
+      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
       !!
       !! history : LIM3.5 - 01-2014 (C. Rousset) original code
       !!-------------------------------------------------------------------
       INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
-
       REAL(wp) ::   zsal, zvi, zvs, zei, zes
       !!-------------------------------------------------------------------
+      at_i (:,:) = 0._wp
+      DO jl = 1, jpl
+         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+      END DO
 
       DO jl = 1, jpl
@@ -534,4 +537,5 @@
                DO ji = 1 , jpi
                   rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
+                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
                   zei              = e_i(ji,jj,jk,jl)
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
@@ -546,4 +550,5 @@
             DO ji = 1 , jpi
                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
+               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
                
                zsal = smv_i(ji,jj,  jl)
@@ -579,8 +584,20 @@
             END DO
          END DO
-      END DO ! jl 
+      END DO 
 
       ! to be sure that at_i is the sum of a_i(jl)
-      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+      at_i (:,:) = 0._wp
+      DO jl = 1, jpl
+         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+      END DO
+
+      ! open water = 1 if at_i=0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
+            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
+         END DO
+      END DO
+
       !
    END SUBROUTINE lim_var_zapsmall
