Index: NEMO/trunk/src/ICE/iceitd.F90
===================================================================
--- NEMO/trunk/src/ICE/iceitd.F90	(revision 10069)
+++ 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
 
