Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5134)
@@ -19,7 +19,5 @@
    USE ice              ! LIM variables
    USE dom_ice          ! LIM domain
-   USE limthd_lac       ! LIM
    USE limvar           ! LIM
-   USE in_out_manager   ! I/O manager
    USE lbclnk           ! lateral boundary condition - MPP exchanges
    USE lib_mpp          ! MPP library
@@ -27,4 +25,5 @@
    USE prtctl           ! Print control
 
+   USE in_out_manager   ! I/O manager
    USE iom              ! I/O manager
    USE lib_fortran      ! glob_sum
@@ -158,5 +157,5 @@
       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
       !-----------------------------------------------------------------------------!
-      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0                ! proport const for PE
+      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0             ! proport const for PE
       !
       CALL lim_itd_me_ridgeprep                                    ! prepare ridging
@@ -275,6 +274,9 @@
          ! 3.4 Compute total area of ice plus open water after ridging.
          !-----------------------------------------------------------------------------!
-
-         CALL lim_itd_me_asumr
+         ! This is in general not equal to one because of divergence during transport
+         asum(:,:) = ato_i(:,:)
+         DO jl = 1, jpl
+            asum(:,:) = asum(:,:) + a_i(:,:,jl)
+         END DO
 
          ! 3.5 Do we keep on iterating ???
@@ -341,19 +343,21 @@
 
       ! Check if there is a ridging error
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
-               WRITE(numout,*) ' '
-               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
-               WRITE(numout,*) ' limitd_me '
-               WRITE(numout,*) ' POINT : ', ji, jj
-               WRITE(numout,*) ' jpl, a_i, athorn '
-               WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
-               DO jl = 1, jpl
-                  WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
-               END DO
-            ENDIF
-         END DO
-      END DO
+      IF( lwp ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
+                  WRITE(numout,*) ' '
+                  WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
+                  WRITE(numout,*) ' limitd_me '
+                  WRITE(numout,*) ' POINT : ', ji, jj
+                  WRITE(numout,*) ' jpl, a_i, athorn '
+                  WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
+                  DO jl = 1, jpl
+                     WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
+                  END DO
+               ENDIF
+            END DO
+         END DO
+      END IF
 
       ! Conservation check
@@ -364,13 +368,13 @@
       ENDIF
 
-      !-----------------------------------------------------------------------------!
-      ! 6) Updating state variables and trend terms (done in limupdate)
-      !-----------------------------------------------------------------------------!
+      ! updates
       CALL lim_var_glo2eqv
       CALL lim_var_zapsmall
       CALL lim_var_agg( 1 ) 
 
-
-      IF(ln_ctl) THEN     ! Control print
+      !-----------------------------------------------------------------------------!
+      ! control prints
+      !-----------------------------------------------------------------------------!
+      IF(ln_ctl) THEN 
          CALL prt_ctl_info(' ')
          CALL prt_ctl_info(' - Cell values : ')
@@ -473,5 +477,5 @@
                      ! PE gain from ridging ice
                      !----------------------------
-                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl)     &
+                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl)     &
                         * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 +  hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )  
                         !!(a**3-b**3)/(a-b) = a*a+ab+b*b                      
@@ -528,12 +532,8 @@
             DO ji = 2, jpim1
                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present
-                  zworka(ji,jj) = 4.0 * strength(ji,jj)              &
-                     &          + strength(ji-1,jj) * tmask(ji-1,jj,1) &  
-                     &          + strength(ji+1,jj) * tmask(ji+1,jj,1) &  
-                     &          + strength(ji,jj-1) * tmask(ji,jj-1,1) &  
-                     &          + strength(ji,jj+1) * tmask(ji,jj+1,1)    
-
-                  zworka(ji,jj) = zworka(ji,jj) /  &
-                     &           ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
+                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
+                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &  
+                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
+                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
                ELSE
                   zworka(ji,jj) = 0._wp
@@ -625,7 +625,9 @@
 
       ! Compute total area of ice plus open water.
-      CALL lim_itd_me_asumr
-      ! This is in general not equal to one 
-      ! because of divergence during transport
+      ! This is in general not equal to one because of divergence during transport
+      asum(:,:) = ato_i(:,:)
+      DO jl = 1, jpl
+         asum(:,:) = asum(:,:) + a_i(:,:,jl)
+      END DO
 
       ! Compute cumulative thickness distribution function
@@ -678,9 +680,9 @@
                DO ji = 1, jpi
                   IF( Gsum(ji,jj,jl) < rn_gstar) THEN
-                     athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
-                        (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
+                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
+                        &                        ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
                   ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
-                     athorn(ji,jj,jl) = Gstari * (rn_gstar-Gsum(ji,jj,jl-1)) *  &
-                        (2.0 - (Gsum(ji,jj,jl-1)+rn_gstar)*Gstari)
+                     athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) *  &
+                        &                        ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
                   ELSE
                      athorn(ji,jj,jl) = 0.0
@@ -693,13 +695,12 @@
          !                        
          zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
-
          DO jl = -1, jpl
             Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
-         END DO !jl
+         END DO
          DO jl = 0, jpl
              athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
          END DO
          !
-      ENDIF ! nn_partfun
+      ENDIF
 
       IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions
@@ -729,5 +730,5 @@
       IF( ln_rafting ) THEN
 
-         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN
+         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN
             DO jl = 1, jpl
                DO jj = 1, jpj
@@ -793,7 +794,7 @@
                ENDIF
 
-            END DO ! ji
-         END DO ! jj
-      END DO ! jl
+            END DO
+         END DO
+      END DO
 
       ! Normalization factor : aksum, ensures mass conservation
@@ -834,5 +835,4 @@
       INTEGER ::   icells            ! number of cells with aicen > puny
       REAL(wp) ::   hL, hR, farea, ztmelts    ! left and right limits of integration
-      REAL(wp) ::   zsstK            ! SST in Kelvin
 
       INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
@@ -910,9 +910,9 @@
                ato_i(ji,jj) = 0._wp
             ENDIF
-         END DO !jj
-      END DO !ji
+         END DO
+      END DO
 
       ! if negative open water area alert it
-      IF( neg_ato_i ) THEN       ! there is a bug
+      IF( neg_ato_i .AND. lwp ) THEN       ! there is a bug
          DO jj = 1, jpj 
             DO ji = 1, jpi
@@ -921,5 +921,5 @@
                   WRITE(numout,*) 'Ridging error: ato_i < 0'
                   WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
-               ENDIF               ! ato_i < -epsi10
+               ENDIF
             END DO
          END DO
@@ -937,5 +937,5 @@
          smv_i_init(:,:,jl) = smv_i(:,:,jl)
          oa_i_init (:,:,jl) = oa_i (:,:,jl)
-      END DO !jl
+      END DO
 
       esnwn_init(:,:,:) = e_s(:,:,1,:)
@@ -968,7 +968,7 @@
                   indxi(icells) = ji
                   indxj(icells) = jj
-               ENDIF ! test on a_icen_init 
-            END DO ! ji
-         END DO ! jj
+               ENDIF
+            END DO
+         END DO
 
          large_afrac = .false.
@@ -1094,5 +1094,5 @@
             dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
 
-         END DO                 ! ij
+         END DO
 
          !--------------------------------------------------------------------
@@ -1112,6 +1112,5 @@
                ! enthalpy of the trapped seawater (J/m2, >0)
                ! clem: if sst>0, then ersw <0 (is that possible?)
-               zsstK  = sst_m(ji,jj) + rt0
-               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) * r1_nlay_i
+               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i
 
                ! heat flux to the ocean
@@ -1121,6 +1120,6 @@
                erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
 
-            END DO ! ij
-         END DO !jk
+            END DO
+         END DO
 
 
@@ -1131,9 +1130,9 @@
                   jj = indxj(ij)
                   eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
-               END DO ! ij
-            END DO !jk
+               END DO
+            END DO
          ENDIF
 
-         IF( large_afrac ) THEN   ! there is a bug
+         IF( large_afrac .AND. lwp ) THEN   ! there is a bug
             DO ij = 1, icells
                ji = indxi(ij)
@@ -1146,5 +1145,5 @@
             END DO
          ENDIF
-         IF( large_afrft ) THEN  ! there is a bug
+         IF( large_afrft .AND. lwp ) THEN  ! there is a bug
             DO ij = 1, icells
                ji = indxi(ij)
@@ -1172,6 +1171,5 @@
                ! Transfer area, volume, and energy accordingly.
 
-               IF( hrmin(ji,jj,jl1) >= hi_max(jl2)  .OR.        &
-                   hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
+               IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
                   hL = 0._wp
                   hR = 0._wp
@@ -1199,5 +1197,5 @@
                   ji = indxi(ij)
                   jj = indxj(ij)
-                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk)
+                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk)
                END DO
             END DO
@@ -1213,6 +1211,5 @@
                ! thickness category jl2, transfer area, volume, and energy accordingly.
                !
-               IF( hraft(ji,jj,jl1) <= hi_max(jl2)  .AND.        &
-                   hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
+               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
                   a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj)
                   v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj)
@@ -1230,6 +1227,5 @@
                   ji = indxi(ij)
                   jj = indxj(ij)
-                  IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)   .AND.        &
-                       hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN
+                  IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1)  ) THEN
                      e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
                   ENDIF
@@ -1271,27 +1267,4 @@
       !
    END SUBROUTINE lim_itd_me_ridgeshift
-
-
-   SUBROUTINE lim_itd_me_asumr
-      !!-----------------------------------------------------------------------------
-      !!                ***  ROUTINE lim_itd_me_asumr ***
-      !!
-      !! ** Purpose :   finds total fractional area
-      !!
-      !! ** Method  :   Find the total area of ice plus open water in each grid cell.
-      !!              This is similar to the aggregate_area subroutine except that the
-      !!              total area can be greater than 1, so the open water area is 
-      !!              included in the sum instead of being computed as a residual. 
-      !!-----------------------------------------------------------------------------
-      INTEGER ::   jl   ! dummy loop index
-      !!-----------------------------------------------------------------------------
-      !
-      asum(:,:) = ato_i(:,:)                    ! open water
-      DO jl = 1, jpl                            ! ice categories
-         asum(:,:) = asum(:,:) + a_i(:,:,jl)
-      END DO
-      !
-   END SUBROUTINE lim_itd_me_asumr
-
 
    SUBROUTINE lim_itd_me_init
@@ -1351,6 +1324,4 @@
    SUBROUTINE lim_itd_me_icestrength
    END SUBROUTINE lim_itd_me_icestrength
-   SUBROUTINE lim_itd_me_sort
-   END SUBROUTINE lim_itd_me_sort
    SUBROUTINE lim_itd_me_init
    END SUBROUTINE lim_itd_me_init
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5134)
@@ -25,5 +25,4 @@
    USE ice              ! LIM-3 variables
    USE limvar           ! LIM-3 variables
-   USE limcons          ! LIM-3 conservation
    USE prtctl           ! Print control
    USE in_out_manager   ! I/O manager
@@ -38,6 +37,4 @@
    PUBLIC   lim_itd_th_rem
    PUBLIC   lim_itd_th_reb
-   PUBLIC   lim_itd_fitline
-   PUBLIC   lim_itd_shiftice
 
    !!----------------------------------------------------------------------
@@ -129,9 +126,10 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes
+               rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes
                ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch
-               rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes
+               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )    !0 if no ice and 1 if yes
                zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch
-               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 
+!clem               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 
+               zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 
             END DO
          END DO
@@ -228,5 +226,5 @@
 
             IF( a_i(ji,jj,kubnd) > epsi10 ) THEN
-               zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1)
+               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) )
             ELSE
                zhbnew(ji,jj,kubnd) = hi_max(kubnd)  
@@ -235,8 +233,6 @@
             ENDIF
 
-            IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1)
-
-         END DO !jj
-      END DO !jj
+         END DO
+      END DO
 
       !-----------------------------------------------------------------------------------------------
@@ -252,10 +248,7 @@
          ij = nind_j(ji) 
 
-         !ji
          IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
             zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
-            ! ji, a_i > epsi10
             IF( zdh0 < 0.0 ) THEN !remove area from category 1
-               ! ji, a_i > epsi10; zdh0 < 0
                zdh0 = MIN( -zdh0, hi_max(klbnd) )
 
@@ -276,8 +269,6 @@
                   v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
                ENDIF
-               ! ji, a_i > epsi10
-
-            ELSE ! if ice accretion
-               ! ji, a_i > epsi10; zdh0 > 0
+
+            ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0
                zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 
                ! zhbnew was 0, and is shifted to the right to account for thin ice
@@ -285,8 +276,7 @@
             ENDIF ! zdh0 
 
-            ! a_i > epsi10
          ENDIF ! a_i > epsi10
 
-      END DO ! ji
+      END DO
 
       !- 7.3 g(h) for each thickness category  
@@ -359,5 +349,5 @@
          ij = nind_j(ji)
          IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN
-            a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
+            a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
             ht_i(ii,ij,1) = rn_himin
          ENDIF
@@ -515,82 +505,83 @@
       END DO
 
-      !----------------------------------------------------------------------------------------------
-      ! 2) Check for daice or dvice out of range, allowing for roundoff error
-      !----------------------------------------------------------------------------------------------
-      ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl
-      ! has a small area, with h(n) very close to a boundary.  Then
-      ! the coefficients of g(h) are large, and the computed daice and
-      ! dvice can be in error. If this happens, it is best to transfer
-      ! either the entire category or nothing at all, depending on which
-      ! side of the boundary hice(n) lies.
-      !-----------------------------------------------------------------
-      DO jl = klbnd, kubnd-1
-
-         zdaice_negative = .false.
-         zdvice_negative = .false.
-         zdaice_greater_aicen = .false.
-         zdvice_greater_vicen = .false.
-
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               IF (zdonor(ji,jj,jl) > 0) THEN
-                  jl1 = zdonor(ji,jj,jl)
-
-                  IF (zdaice(ji,jj,jl) < 0.0) THEN
-                     IF (zdaice(ji,jj,jl) > -epsi10) THEN
-                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  &
-                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
-                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category
-                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
-                        ELSE
-                           zdaice(ji,jj,jl) = 0.0 ! shift no ice
-                           zdvice(ji,jj,jl) = 0.0
-                        ENDIF
-                     ELSE
-                        zdaice_negative = .true.
-                     ENDIF
-                  ENDIF
-
-                  IF (zdvice(ji,jj,jl) < 0.0) THEN
-                     IF (zdvice(ji,jj,jl) > -epsi10 ) THEN
-                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  &
-                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
-                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
-                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
-                        ELSE
-                           zdaice(ji,jj,jl) = 0.0    ! shift no ice
-                           zdvice(ji,jj,jl) = 0.0
-                        ENDIF
-                     ELSE
-                        zdvice_negative = .true.
-                     ENDIF
-                  ENDIF
-
-                  ! If daice is close to aicen, set daice = aicen.
-                  IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN
-                     IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN
-                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
-                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
-                     ELSE
-                        zdaice_greater_aicen = .true.
-                     ENDIF
-                  ENDIF
-
-                  IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN
-                     IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN
-                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
-                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
-                     ELSE
-                        zdvice_greater_vicen = .true.
-                     ENDIF
-                  ENDIF
-
-               ENDIF               ! donor > 0
-            END DO
-         END DO
-
-      END DO
-
+!clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10)
+!      !----------------------------------------------------------------------------------------------
+!      ! 2) Check for daice or dvice out of range, allowing for roundoff error
+!      !----------------------------------------------------------------------------------------------
+!      ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl
+!      ! has a small area, with h(n) very close to a boundary.  Then
+!      ! the coefficients of g(h) are large, and the computed daice and
+!      ! dvice can be in error. If this happens, it is best to transfer
+!      ! either the entire category or nothing at all, depending on which
+!      ! side of the boundary hice(n) lies.
+!      !-----------------------------------------------------------------
+!      DO jl = klbnd, kubnd-1
+!
+!         zdaice_negative = .false.
+!         zdvice_negative = .false.
+!         zdaice_greater_aicen = .false.
+!         zdvice_greater_vicen = .false.
+!
+!         DO jj = 1, jpj
+!            DO ji = 1, jpi
+!
+!               IF (zdonor(ji,jj,jl) > 0) THEN
+!                  jl1 = zdonor(ji,jj,jl)
+!
+!                  IF (zdaice(ji,jj,jl) < 0.0) THEN
+!                     IF (zdaice(ji,jj,jl) > -epsi10) THEN
+!                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  &
+!                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
+!                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category
+!                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
+!                        ELSE
+!                           zdaice(ji,jj,jl) = 0.0 ! shift no ice
+!                           zdvice(ji,jj,jl) = 0.0
+!                        ENDIF
+!                     ELSE
+!                        zdaice_negative = .true.
+!                     ENDIF
+!                  ENDIF
+!
+!                  IF (zdvice(ji,jj,jl) < 0.0) THEN
+!                     IF (zdvice(ji,jj,jl) > -epsi10 ) THEN
+!                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  &
+!                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
+!                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
+!                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
+!                        ELSE
+!                           zdaice(ji,jj,jl) = 0.0    ! shift no ice
+!                           zdvice(ji,jj,jl) = 0.0
+!                        ENDIF
+!                    ELSE
+!                       zdvice_negative = .true.
+!                    ENDIF
+!                 ENDIF
+!
+!                  ! If daice is close to aicen, set daice = aicen.
+!                  IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN
+!                     IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN
+!                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
+!                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
+!                    ELSE
+!                       zdaice_greater_aicen = .true.
+!                    ENDIF
+!                  ENDIF
+!
+!                  IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN
+!                     IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN
+!                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
+!                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
+!                     ELSE
+!                        zdvice_greater_vicen = .true.
+!                     ENDIF
+!                  ENDIF
+!
+!               ENDIF               ! donor > 0
+!            END DO
+!         END DO
+!
+!      END DO
+!clem
       !-------------------------------------------------------------------------------
       ! 3) Transfer volume and energy between categories
@@ -614,6 +605,6 @@
 
             jl1 = zdonor(ii,ij,jl)
-            rswitch       = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) )
-            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch
+            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi20 ) )
+            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi20 ) * rswitch
             IF( jl1 == jl) THEN   ;   jl2 = jl1+1
             ELSE                  ;   jl2 = jl 
@@ -710,5 +701,4 @@
                   ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl) 
                   t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 
-                  rswitch         =  1.0 - MAX( 0.0, SIGN( 1.0, -v_s(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes
                ELSE
                   ht_i(ji,jj,jl)  = 0._wp
@@ -765,9 +755,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi 
-               IF( a_i(ji,jj,jl) > epsi10 ) THEN 
-                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-               ELSE
-                  ht_i(ji,jj,jl) = 0._wp
-               ENDIF
+               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
+               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
             END DO
          END DO
@@ -775,19 +762,5 @@
 
       !------------------------------------------------------------------------------
-      ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd)
-      !------------------------------------------------------------------------------
-      DO jj = 1, jpj 
-         DO ji = 1, jpi 
-            IF( a_i(ji,jj,klbnd) > epsi10 ) THEN
-               IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN
-                  a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max(0) 
-                  ht_i(ji,jj,klbnd) = hi_max(0)
-               ENDIF
-            ENDIF
-         END DO
-      END DO
-
-      !------------------------------------------------------------------------------
-      ! 3) If a category thickness is not in bounds, shift the
+      ! 2) If a category thickness is not in bounds, shift the
       ! entire area, volume, and energy to the neighboring category
       !------------------------------------------------------------------------------
@@ -818,12 +791,10 @@
                   zdonor(ji,jj,jl)  = jl 
                   ! begin TECLIM change
-                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl)
-                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)
                   !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp
                   !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp
                   ! end TECLIM change 
                   ! clem: how much of a_i you send in cat sup is somewhat arbitrary
-                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi10 ) / ht_i(ji,jj,jl)  
-                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 )
+                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)  
+                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 )
                ENDIF
             END DO
@@ -839,5 +810,5 @@
          ENDIF
          !
-      END DO                    ! jl
+      END DO
 
       !----------------------------
@@ -888,5 +859,5 @@
 
       !------------------------------------------------------------------------------
-      ! 4) Conservation check
+      ! 3) Conservation check
       !------------------------------------------------------------------------------
 
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5134)
@@ -115,5 +115,5 @@
                DO ji = 1, jpi
                   !0 if no ice and 1 if yes
-                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 )  )
+                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )
                   !Energy of melting q(S,T) [J.m-3]
                   e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
@@ -125,5 +125,5 @@
                DO ji = 1, jpi
                   !0 if no ice and 1 if yes
-                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi20 )  )
+                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 )  )
                   !Energy of melting q(S,T) [J.m-3]
                   e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
@@ -158,5 +158,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            rswitch  = tmask(ji,jj,1) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice
+            rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
             !
             !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
@@ -360,4 +360,5 @@
       CALL lim_var_eqv2glo
 
+      CALL lim_var_zapsmall
       !--------------------------------------------
       ! Diagnostic thermodynamic growth rates
@@ -413,8 +414,8 @@
 
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !------------------------------------------------------------------------------|
       !  7) Add frazil ice growing in leads.
       !------------------------------------------------------------------------------|
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       CALL lim_thd_lac
       CALL lim_var_glo2eqv    ! only for info
@@ -513,5 +514,5 @@
             a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 
             ! adjust thickness
-            rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) )
+            rswitch     = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) )
             ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 )
             ! retrieve total concentration
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5134)
@@ -244,5 +244,5 @@
          ! thickness change
          IF( zdh_s_pre(ji) > 0._wp ) THEN
-         rswitch        = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )
+         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
          zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
          zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 
@@ -266,5 +266,5 @@
             ! thickness change
             rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
-            rswitch          = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) ) 
+            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 
             zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting
@@ -321,7 +321,7 @@
       DO jk = 1, nlay_s
          DO ji = kideb,kiut
-            rswitch       =  MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 )  )
-            q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) *             &
-              &            ( (   MAX( 0._wp, dh_s_tot(ji) )              ) * zqprec(ji) +  &
+            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  )
+            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          &
+              &            ( (   MAX( 0._wp, dh_s_tot(ji) )               ) * zqprec(ji) +  &
               &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
             zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
@@ -453,7 +453,7 @@
                q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
                
-            ENDIF ! fc_bo_i
-         END DO ! ji
-      END DO ! iter
+            ENDIF
+         END DO
+      END DO
 
       ! Contribution to Energy and Salt Fluxes
@@ -592,6 +592,6 @@
          zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
          rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
-         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, zq_s(ji) - epsi20 ) )
-         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( zq_s(ji), epsi20 )
+         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )
+         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )
          zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
          zdh_s_mel(ji)   = zdh_s_mel(ji) + zdeltah(ji,1)    
@@ -599,7 +599,7 @@
          ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1)
         
-         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji)                ! update available heat (J.m-2)
+         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2)
          ! heat used to melt snow
-         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
+         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
          ! Contribution to mass flux
          wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
@@ -630,12 +630,12 @@
 
          ! entrapment during snow ice formation
-         ! new salinity difference stored (to be used in limthd_ent.F90)
+         ! new salinity difference stored (to be used in limthd_sal.F90)
          IF (  nn_icesal == 2  ) THEN
-            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) )
+            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
             ! salinity dif due to snow-ice formation
-            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch     
+            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch     
             ! salinity dif due to bottom growth 
             IF (  zf_tt(ji)  < 0._wp ) THEN
-               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch
+               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch
             ENDIF
          ENDIF
@@ -663,8 +663,5 @@
          h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
          
-         ! Total ablation (to debug)
-         IF( ht_i_1d(ji) <= 0._wp )   a_i_1d(ji) = 0._wp
-
-      END DO !ji
+      END DO
 
       !
@@ -676,15 +673,18 @@
          rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
          t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
-      END DO  ! ji
+      END DO
 
       DO jk = 1, nlay_s
          DO ji = kideb,kiut
             ! mask enthalpy
-            rswitch       =  MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
-            q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk)
+            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
+            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)
             ! recalculate t_s_1d from q_s_1d
-            t_s_1d(ji,jk) = rt0 + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
-         END DO
-      END DO
+            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
+         END DO
+      END DO
+
+      ! --- ensure that a_i = 0 where ht_i = 0 ---
+      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
       
       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema )
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 5134)
@@ -132,5 +132,5 @@
       DO jk1 = 1, nlay_i
          DO ji = kideb, kiut
-            rswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi20 ) ) 
+            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 
             qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
          ENDDO
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5134)
@@ -130,5 +130,5 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi20 )  )   !0 if no ice
+                  rswitch          = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )   !0 if no ice
                   e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp )
                END DO
@@ -482,5 +482,5 @@
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               rswitch          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes
+               rswitch          = MAX( 0._wp , SIGN( 1._wp , za_i_1d(ji,jl) - epsi20 ) )  ! 0 if no ice and 1 if yes
                zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch   
             END DO 
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5134)
@@ -68,8 +68,8 @@
       CHARACTER(len=80) ::   cltmp
       !
-      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm, zs0at
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ow
-      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
+      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
+      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold, zsmvold  ! old ice volume...
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
@@ -80,8 +80,8 @@
       IF( nn_timing == 1 )  CALL timing_start('limtrp')
 
-      CALL wrk_alloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold )
-      CALL wrk_alloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e )
-      CALL wrk_alloc( jpi,jpj,1,         zs0ow )
-      CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, zs0e )
+      CALL wrk_alloc( jpi,jpj,           zsm, zatold, zeiold, zesold )
+      CALL wrk_alloc( jpi,jpj,jpl,       z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+      CALL wrk_alloc( jpi,jpj,1,         z0opw )
+      CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, z0ei )
       CALL wrk_alloc( jpi,jpj,jpl,       zhimax, zviold, zvsold, zsmvold )
 
@@ -118,9 +118,9 @@
          ! in case advection creates ice too thick.
          !---------------------------------------------------------------------
-         zhimax(:,:,:) = ht_i(:,:,:)
+         zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
          DO jl = 1, jpl
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
+                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
                END DO
             END DO
@@ -155,14 +155,14 @@
          ! transported fields                                        
          !-------------------------
-         zs0ow(:,:,1) = ato_i(:,:) * e12t(:,:)              ! Open water area 
-         DO jl = 1, jpl
-            zs0sn (:,:,jl)   = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
-            zs0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
-            zs0a  (:,:,jl)   = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
-            zs0sm (:,:,jl)   = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
-            zs0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
-            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
+         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area 
+         DO jl = 1, jpl
+            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
+            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
+            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
+            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
+            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
+            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
             DO jk = 1, nlay_i
-               zs0e  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
+               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
             END DO
          END DO
@@ -171,38 +171,38 @@
          IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
             DO jt = 1, initad
-               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area
+               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
                   &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &
+               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
                   &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
                DO jl = 1, jpl
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
                   DO jk = 1, nlay_i                                                                !--- ice heat contents ---
-                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
@@ -212,39 +212,39 @@
          ELSE
             DO jt = 1, initad
-               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area
+               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
                   &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &
+               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
                   &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
                DO jl = 1, jpl
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
 
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
                   DO jk = 1, nlay_i                                                           !--- ice heat contents ---
-                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
@@ -257,14 +257,14 @@
          ! Recover the properties from their contents
          !-------------------------------------------
-         ato_i(:,:) = zs0ow(:,:,1) * r1_e12t(:,:)
-         DO jl = 1, jpl
-            v_i  (:,:,jl)   = zs0ice(:,:,jl) * r1_e12t(:,:)
-            v_s  (:,:,jl)   = zs0sn (:,:,jl) * r1_e12t(:,:)
-            smv_i(:,:,jl)   = zs0sm (:,:,jl) * r1_e12t(:,:)
-            oa_i (:,:,jl)   = zs0oi (:,:,jl) * r1_e12t(:,:)
-            a_i  (:,:,jl)   = zs0a  (:,:,jl) * r1_e12t(:,:)
-            e_s  (:,:,1,jl) = zs0c0 (:,:,jl) * r1_e12t(:,:)
+         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
+         DO jl = 1, jpl
+            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
+            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
+            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
+            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
+            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
+            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
             DO jk = 1, nlay_i
-               e_i  (:,:,jk,jl) = zs0e  (:,:,jk,jl) * r1_e12t(:,:)
+               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
             END DO
          END DO
@@ -293,5 +293,5 @@
          END DO
          !
-         CALL lim_hdf( ato_i (:,:) )   ! Diffusion
+         CALL lim_hdf( ato_i (:,:) )
 
          !------------------------------------
@@ -356,4 +356,5 @@
 
                   IF ( v_i(ji,jj,jl) > 0._wp ) THEN
+
                      zvi  = v_i  (ji,jj,jl)
                      zvs  = v_s  (ji,jj,jl)
@@ -361,31 +362,29 @@
                      zes  = e_s  (ji,jj,1,jl)
                      zei  = SUM( e_i(ji,jj,1:nlay_i,jl) )
-                     zdv  = v_i(ji,jj,jl) - zviold(ji,jj,jl)   
-                     
-                     rswitch = 1._wp
-                     IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
-                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                          
-                        ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
-                        rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
-                        a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
-                     ELSE
-                        ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
-                        rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
-                        a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
+
+                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)  
+
+                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
+                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN                                          
+
+                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
+                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
+
+                        ! small correction due to *rswitch for a_i
+                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
+                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
+                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
+                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
+                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
+
+                        ! Update mass fluxes
+                        wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
+                        wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
+                        sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
+                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
+                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
+
                      ENDIF
 
-                     ! small correction due to *rswitch for a_i
-                     v_i  (ji,jj,jl) = rswitch * v_i  (ji,jj,jl)
-                     v_s  (ji,jj,jl) = rswitch * v_s  (ji,jj,jl)
-                     smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)
-                     e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)
-                     e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
-
-                     ! Update mass fluxes
-                     wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
-                     wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
-                     sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
-                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
-                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
                   ENDIF
 
@@ -423,21 +422,17 @@
          vt_s (:,:) = 0._wp
          at_i (:,:) = 0._wp
-         !
          DO jl = 1, jpl
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  !
-                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
-                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
-                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
-               END DO
-            END DO
-         END DO
-         ! -------------------------------------------------
-
-         ! open water
+                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
+                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
+                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
+               END DO
+            END DO
+         END DO
+
+         ! --- open water = 1 if at_i=0 --------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ! open water = 1 if at_i=0
                rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
                ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
@@ -457,8 +452,8 @@
       IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
       !
-      CALL wrk_dealloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold )
-      CALL wrk_dealloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e )
-      CALL wrk_dealloc( jpi,jpj,1,         zs0ow )
-      CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, zs0e )
+      CALL wrk_dealloc( jpi,jpj,           zsm, zatold, zeiold, zesold )
+      CALL wrk_dealloc( jpi,jpj,jpl,       z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+      CALL wrk_dealloc( jpi,jpj,1,         z0opw )
+      CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, z0ei )
       CALL wrk_dealloc( jpi,jpj,jpl,       zviold, zvsold, zhimax, zsmvold )
       !
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5134)
@@ -69,16 +69,8 @@
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      !-----------------
-      ! zap small values
-      !-----------------
-      CALL lim_var_zapsmall
-
       CALL lim_var_glo2eqv
-     
       !----------------------------------------------------
-      ! Rebin categories with thickness out of bounds
-      !----------------------------------------------------
-      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl)
-
+      ! ice concentration should not exceed amax 
+      !-----------------------------------------------------
       at_i(:,:) = 0._wp
       DO jl = 1, jpl
@@ -86,7 +78,4 @@
       END DO
 
-      !----------------------------------------------------
-      ! ice concentration should not exceed amax 
-      !-----------------------------------------------------
       DO jl  = 1, jpl
          DO jj = 1, jpj
@@ -94,18 +83,12 @@
                IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
                   a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
-                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
                ENDIF
             END DO
          END DO
       END DO
-
-      at_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
     
-      ! --------------------------------------
-      ! Final thickness distribution rebinning
-      ! --------------------------------------
+      !----------------------------------------------------
+      ! Rebin categories with thickness out of bounds
+      !----------------------------------------------------
       IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)
 
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5134)
@@ -56,5 +56,5 @@
       INTEGER, INTENT(in) ::   kt    ! number of iteration
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zh, zsal
+      REAL(wp) ::   zsal
       REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-------------------------------------------------------------------
@@ -69,26 +69,12 @@
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      !-----------------
-      ! zap small values
-      !-----------------
-      CALL lim_var_agg( 1 )
-      CALL lim_var_zapsmall
-      CALL lim_var_glo2eqv
-
-      !----------------------------------------------------
-      ! Rebin categories with thickness out of bounds
-      !----------------------------------------------------
-      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl)
-
       !----------------------------------------------------------------------
       ! Constrain the thickness of the smallest category above himin
       !----------------------------------------------------------------------
+      CALL lim_var_glo2eqv
       DO jj = 1, jpj 
          DO ji = 1, jpi
             IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
-               zh             = rn_himin / ht_i(ji,jj,1)
-               ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh
-               ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh
-               a_i (ji,jj,1) = a_i(ji,jj,1)  / zh
+               a_i (ji,jj,1) = a_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin
             ENDIF
          END DO
@@ -108,5 +94,4 @@
                IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
                   a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
-                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
                ENDIF
             END DO
@@ -114,12 +99,7 @@
       END DO
 
-      at_i(:,:) = 0.0
-      DO jl = 1, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-
-      ! --------------------------------------
-      ! Final thickness distribution rebinning
-      ! --------------------------------------
+      !----------------------------------------------------
+      ! Rebin categories with thickness out of bounds
+      !----------------------------------------------------
       IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
 
@@ -140,10 +120,10 @@
                   ! salinity stays in bounds
                   rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
-                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) !+ rn_simin * ( 1._wp - rswitch ) * v_i(ji,jj,jl)
+                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
                   ! associated salt flux
                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
-               END DO ! ji
-            END DO ! jj
-         END DO !jl
+               END DO
+            END DO
+         END DO
       ENDIF
 
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5133)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5134)
@@ -161,5 +161,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
+               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )   !0 if no ice and 1 if yes
                ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
@@ -173,5 +173,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
+                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) )   !0 if no ice and 1 if yes
                   sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch
                END DO
@@ -190,5 +190,5 @@
                DO ji = 1, jpi
                   !                                                              ! Energy of melting q(S,T) [J.m-3]
-                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
+                  rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
                   zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 
                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0              ! Ice layer melt temperature
@@ -215,5 +215,5 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
+                  rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
                   zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
                   !
@@ -233,5 +233,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
+                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
                   tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
@@ -339,8 +339,8 @@
                      !                                      ! weighting the profile
                      s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
-                  END DO ! ji
-               END DO ! jj
-            END DO ! jk
-         END DO ! jl
+                  END DO
+               END DO
+            END DO
+         END DO
          !
       ENDIF ! nn_icesal
@@ -384,5 +384,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
+                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
                   tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
                      &                      * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 )
@@ -572,8 +572,4 @@
                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
 
-               ! ice salinity must stay in bounds
-               IF(  nn_icesal == 2  ) THEN
-                  smv_i(ji,jj,jl) = MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
-               ENDIF
                ! update exchanges with ocean
                sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
