Index: /NEMO/trunk/src/ICE/icecor.F90
===================================================================
--- /NEMO/trunk/src/ICE/icecor.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icecor.F90	(revision 10994)
@@ -66,7 +66,6 @@
          WRITE(numout,*) '~~~~~~~'
       ENDIF
-      !
       !                             !-----------------------------------------------------
-      !                             !  ice thickness must exceed himin (for ice diff)    !
+      !                             !  ice thickness must exceed himin (for temp. diff.) !
       !                             !-----------------------------------------------------
       WHERE( a_i(:,:,:) >= epsi20 )   ;   h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
@@ -97,5 +96,4 @@
          END DO
       ENDIF
-
       !                             !-----------------------------------------------------
       !                             !  Rebin categories with thickness out of bounds     !
Index: /NEMO/trunk/src/ICE/icectl.F90
===================================================================
--- /NEMO/trunk/src/ICE/icectl.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icectl.F90	(revision 10994)
@@ -69,5 +69,5 @@
       !!
       REAL(wp) ::   zv, zs, zt, zfs, zfv, zft
-      REAL(wp) ::   zvmin, zamin, zamax 
+      REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin
       REAL(wp) ::   zvtrp, zetrp
       REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill
@@ -141,7 +141,10 @@
          zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv
 
-         zvmin = glob_min( 'icectl', v_i )
-         zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) )
-         zamin = glob_min( 'icectl', a_i )
+         zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) )
+         zvmin  = glob_min( 'icectl', v_i )
+         zamin  = glob_min( 'icectl', a_i )
+         zsmin  = glob_min( 'icectl', sv_i )
+         zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) )
+         zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) )
 
          ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
@@ -152,13 +155,18 @@
 
          IF(lwp) THEN
-            IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv
-            IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs
-            IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt
-            IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
-            IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10   &
-               & .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' .AND. cd_routine /= 'Hbig' ) &
-               &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
-            IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
-!clem: the following check fails when using UMx advection scheme (see comments in icedyn_adv.F90)
+            ! check conservation issues
+            IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv
+            IF ( ABS( zs ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs
+            IF ( ABS( zt ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt
+            ! check maximum ice concentration
+            IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  &
+               &                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
+            ! check negative values
+            IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
+            IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+            IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin
+            IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin
+            IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin
+!clem: the following check fails (I think...)
 !            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN
 !                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
Index: /NEMO/trunk/src/ICE/icedyn.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icedyn.F90	(revision 10994)
@@ -118,4 +118,5 @@
          CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
          CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
+         CALL ice_var_zapsmall                                              ! -- zap small areas
          !
       CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
@@ -189,6 +190,4 @@
       ! controls
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-      CALL ice_var_zapsmall                       !-- zap small areas
       !
       at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
Index: /NEMO/trunk/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10994)
@@ -142,8 +142,4 @@
       INTEGER, PARAMETER ::   jp_itermax = 20    
       !!-------------------------------------------------------------------
-      ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem)
-      !       likely due to truncation error ( i.e. 1. - 1. /= 0 )
-      !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
-      
       ! controls
       IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
@@ -156,14 +152,14 @@
       ENDIF      
 
-      CALL ice_var_zapsmall   ! Zero out categories with very small areas
-
       !--------------------------------
       ! 0) Identify grid cells with ice
       !--------------------------------
+      at_i(:,:) = SUM( a_i, dim=3 )
+      !
       npti = 0   ;   nptidx(:) = 0
       ipti = 0   ;   iptidx(:) = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF ( at_i(ji,jj) > 0._wp ) THEN
+            IF ( at_i(ji,jj) > epsi10 ) THEN
                npti           = npti + 1
                nptidx( npti ) = (jj - 1) * jpi + ji
@@ -178,10 +174,10 @@
         
          ! just needed here
-         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i  )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i )
          ! needed here and in the iteration loop
-         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) )
-         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i )
 
          DO ji = 1, npti
@@ -310,6 +306,6 @@
 
       !                       ! Ice thickness needed for rafting
-      WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
-      ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp
+      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
       END WHERE
 
@@ -329,6 +325,6 @@
       zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
       !
-      WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
-      ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp
+      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
+      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
       END WHERE
       !
@@ -455,6 +451,6 @@
       ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.  
       ! NOTE: 0 < aksum <= 1
-      WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
-      ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp
+      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
+      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
       END WHERE
       
@@ -466,5 +462,5 @@
          DO ji = 1, npti
             zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
-            IF( zfac > pa_i(ji,jl) ) THEN
+            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
                closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
             ENDIF
@@ -510,4 +506,5 @@
       REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
       REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
+      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i)
       !
       REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
@@ -518,5 +515,7 @@
       INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
       !!-------------------------------------------------------------------
-
+      !
+      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume
+      !
       ! 1) Change in open water area due to closing and opening
       !--------------------------------------------------------
@@ -535,6 +534,8 @@
             IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
 
-               z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
-
+               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
+               ELSE                                 ;   z1_ai(ji) = 0._wp
+               ENDIF
+               
                ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
                airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
@@ -549,10 +550,13 @@
 
                ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
-               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
+               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg
+               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0
+               ELSE                           ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0
+               ENDIF
                ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
 
                ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
                virdg1     = v_i_2d (ji,jl1)   * afrdg
-               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg )
+               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw
                vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
                sirdg1     = sv_i_2d(ji,jl1)   * afrdg
@@ -726,7 +730,8 @@
       END DO ! jl1
       !
+      ! roundoff errors
+      !----------------
       ! In case ridging/rafting lead to very small negative values (sometimes it happens)
-      WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp
-      WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
       !
    END SUBROUTINE rdgrft_shift
@@ -854,5 +859,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
-
+         !
          !                 !---------------------!
       CASE( 2 )            !==  from 1D to 2D  ==!
Index: /NEMO/trunk/src/ICE/iceitd.F90
===================================================================
--- /NEMO/trunk/src/ICE/iceitd.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/iceitd.F90	(revision 10994)
@@ -21,4 +21,5 @@
    USE ice1D          ! sea-ice: thermodynamic variables
    USE ice            ! sea-ice: variables
+   USE icevar         ! sea-ice: operations
    USE icectl         ! sea-ice: conservation tests
    USE icetab         ! sea-ice: convert 1D<=>2D
@@ -91,4 +92,6 @@
       !  1) Identify grid cells with ice
       !-----------------------------------------------------------------------------------------------
+      at_i(:,:) = SUM( a_i, dim=3 )
+      !
       npti = 0   ;   nptidx(:) = 0
       DO jj = 1, jpj
@@ -249,5 +252,5 @@
             ! --- g(h) for each thickness category --- !  
             CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
-               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR   (1:npti,jl)   )   ! out
+               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out
             !
          END DO
@@ -389,9 +392,11 @@
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
       !
-      INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices
-      INTEGER  ::   ii, ij, jl2, jl1   ! local integers
+      INTEGER  ::   ji, jl, jk         ! dummy loop indices
+      INTEGER  ::   jl2, jl1           ! local integers
       REAL(wp) ::   ztrans             ! ice/snow transferred
-      REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace
-      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    -
+      REAL(wp), DIMENSION(jpij)            ::   zworka, zworkv   ! workspace
+      REAL(wp), DIMENSION(jpij,jpl)        ::   zaTsfn           !  -    -
+      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d
+      REAL(wp), DIMENSION(jpij,nlay_s,jpl) ::   ze_s_2d
       !!------------------------------------------------------------------
          
@@ -405,4 +410,14 @@
       CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
       CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
+      DO jl = 1, jpl
+         DO jk = 1, nlay_s
+            CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
+         END DO
+         DO jk = 1, nlay_i
+            CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+         END DO
+      END DO
+      ! to correct roundoff errors on a_i
+      CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d )
 
       !----------------------------------------------------------------------------------------------
@@ -435,8 +450,4 @@
                ELSE                                  ;   zworka(ji) = 0._wp
                ENDIF
-               !
-               ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
-               !       because of truncation error ( i.e. 1. - 1. /= 0 )
-               !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
                !
                a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
@@ -476,8 +487,5 @@
          !
          DO jk = 1, nlay_s         !--- Snow heat content
-            !
             DO ji = 1, npti
-               ii = MOD( nptidx(ji) - 1, jpi ) + 1
-               ij = ( nptidx(ji) - 1 ) / jpi + 1
                !
                jl1 = kdonor(ji,jl)
@@ -487,8 +495,7 @@
                   ELSE                ;  jl2 = jl
                   ENDIF
-                  !
-                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
-                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
-                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
+                  ztrans             = ze_s_2d(ji,jk,jl1) * zworkv(ji)
+                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans
+                  ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans
                ENDIF
             END DO
@@ -497,6 +504,4 @@
          DO jk = 1, nlay_i         !--- Ice heat content
             DO ji = 1, npti
-               ii = MOD( nptidx(ji) - 1, jpi ) + 1
-               ij = ( nptidx(ji) - 1 ) / jpi + 1
                !
                jl1 = kdonor(ji,jl)
@@ -506,8 +511,7 @@
                   ELSE                ;  jl2 = jl
                   ENDIF
-                  !
-                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
-                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
-                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
+                  ztrans             = ze_i_2d(ji,jk,jl1) * zworkv(ji)
+                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans
+                  ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans
                ENDIF
             END DO
@@ -515,7 +519,21 @@
          !
       END DO                   ! boundaries, 1 to jpl-1
+
+      !-------------------
+      ! 3) roundoff errors
+      !-------------------
+      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
+      !       because of truncation error ( i.e. 1. - 1. /= 0 )
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
+
+      ! at_i must be <= rn_amax
+      zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
+      DO jl  = 1, jpl
+         WHERE( zworka(1:npti) > rn_amax_1d(1:npti) )   &
+            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti)
+      END DO
       
       !-------------------------------------------------------------------------------
-      ! 3) Update ice thickness and temperature
+      ! 4) Update ice thickness and temperature
       !-------------------------------------------------------------------------------
       WHERE( a_i_2d(1:npti,:) >= epsi20 )
@@ -536,4 +554,12 @@
       CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
       CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
+      DO jl = 1, jpl
+         DO jk = 1, nlay_s
+            CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
+         END DO
+         DO jk = 1, nlay_i
+            CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+         END DO
+      END DO
       !
    END SUBROUTINE itd_shiftice
@@ -558,4 +584,6 @@
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
+      !
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
       !
       jdonor(:,:) = 0
@@ -635,4 +663,6 @@
       END DO
       !
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      !
    END SUBROUTINE ice_itd_reb
 
Index: /NEMO/trunk/src/ICE/icestp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icestp.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icestp.F90	(revision 10994)
@@ -323,4 +323,7 @@
          WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s
       ENDIF
+      !                                        !--- change max ice concentration for roundoff errors
+      rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 )
+      rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 )
       !                                        !--- check consistency
       IF ( jpl > 1 .AND. ln_virtual_itd ) THEN
Index: /NEMO/trunk/src/ICE/icethd.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icethd.F90	(revision 10994)
@@ -102,6 +102,4 @@
       ENDIF
       
-      CALL ice_var_glo2eqv
-
       !---------------------------------------------!
       ! computation of friction velocity at T points
@@ -162,6 +160,7 @@
             qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(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 
-            IF( zqld > 0._wp ) THEN
+            ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 
+            ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget
+            IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN
                fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
                qlead(ji,jj) = 0._wp
@@ -242,9 +241,6 @@
          !
       END DO
-
+      !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-      !
-                              CALL ice_var_zapsmall                 ! --- remove very small ice concentration (<1e-10) --- !
-      !                                                             !     & make sure at_i=sum(a_i) & ato_i=1 where at_i=0
       !                   
       IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
Index: /NEMO/trunk/src/ICE/icethd_do.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_do.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icethd_do.F90	(revision 10994)
@@ -114,7 +114,5 @@
       IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
 
-      CALL ice_var_agg(1)
-      CALL ice_var_glo2eqv
-
+      at_i(:,:) = SUM( a_i, dim=3 )
       !------------------------------------------------------------------------------!
       ! 1) Collection thickness of ice formed in leads and polynyas
@@ -319,12 +317,12 @@
 
          ! --- lateral ice growth --- !
-         ! If lateral ice growth gives an ice concentration gt 1, then
+         ! If lateral ice growth gives an ice concentration > amax, then
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, npti
-            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) )
+            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
+               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
-               za_newice(ji) = za_newice(ji) - zda_res  (ji)
-               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
+               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) )
+               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) )
             ELSE
                zda_res(ji) = 0._wp
Index: /NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- /NEMO/trunk/src/ICE/icevar.F90	(revision 10993)
+++ /NEMO/trunk/src/ICE/icevar.F90	(revision 10994)
@@ -44,5 +44,6 @@
    !!   ice_var_salprof1d : salinity profile in the ice 1D
    !!   ice_var_zapsmall  : remove very small area and volume
-   !!   ice_var_zapneg    : remove negative ice fields (to debug the advection scheme UM3-5)
+   !!   ice_var_zapneg    : remove negative ice fields
+   !!   ice_var_roundoff  : remove negative values arising from roundoff erros
    !!   ice_var_itd       : convert 1-cat to jpl-cat
    !!   ice_var_itd2      : convert N-cat to jpl-cat
@@ -71,4 +72,5 @@
    PUBLIC   ice_var_zapsmall
    PUBLIC   ice_var_zapneg
+   PUBLIC   ice_var_roundoff
    PUBLIC   ice_var_itd
    PUBLIC   ice_var_itd2
@@ -622,4 +624,36 @@
    END SUBROUTINE ice_var_zapneg
 
+
+   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_roundoff ***
+      !!
+      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
+      !!-------------------------------------------------------------------
+      !
+      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
+      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
+      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
+      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
+      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
+      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
+      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
+      IF ( ln_pnd_H12 ) THEN
+         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
+         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
+      ENDIF
+      !
+   END SUBROUTINE ice_var_roundoff
+   
    
    SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
@@ -656,6 +690,6 @@
       INTEGER  :: idim, i_fill, jl0  
       REAL(wp) :: zarg, zV, zconv, zdh, zdv
-      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
+      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
       INTEGER , DIMENSION(4)                  ::   itest
       !!-------------------------------------------------------------------
@@ -786,5 +820,5 @@
       !!
       !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
-      !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
+       !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
       !!              
       !!               3) Expand the filling to the empty cat between jlmin and jlmax 
