Index: NEMO/releases/release-4.0/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/icedyn_rdgrft.F90	(revision 10944)
+++ NEMO/releases/release-4.0/src/ICE/icedyn_rdgrft.F90	(revision 10993)
@@ -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  ==!
