Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 8342)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 8355)
@@ -28,7 +28,4 @@
    PRIVATE
 
-   PUBLIC   lim_column_sum
-   PUBLIC   lim_column_sum_energy
-   PUBLIC   lim_cons_check
    PUBLIC   lim_cons_hsm
    PUBLIC   lim_cons_final
@@ -40,120 +37,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE lim_column_sum( ksum, pin, pout )
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE lim_column_sum ***
-      !!
-      !! ** Purpose : Compute the sum of xin over nsum categories
-      !!
-      !! ** Method  : Arithmetics
-      !!
-      !! ** Action  : Gets xin(ji,jj,jl) and computes xout(ji,jj)
-      !!---------------------------------------------------------------------
-      INTEGER                   , INTENT(in   ) ::   ksum   ! number of categories/layers
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pin    ! input field
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   pout   ! output field
-      !
-      INTEGER ::   jl   ! dummy loop indices
-      !!---------------------------------------------------------------------
-      !
-      pout(:,:) = pin(:,:,1)
-      DO jl = 2, ksum
-         pout(:,:) = pout(:,:) + pin(:,:,jl)
-      END DO
-      !
-   END SUBROUTINE lim_column_sum
-
-
-   SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout)
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE lim_column_sum_energy ***
-      !!
-      !! ** Purpose : Compute the sum of xin over nsum categories
-      !!              and nlay layers
-      !!
-      !! ** Method  : Arithmetics
-      !!---------------------------------------------------------------------
-      INTEGER                                , INTENT(in   ) ::   ksum   !: number of categories
-      INTEGER                                , INTENT(in   ) ::   klay   !: number of vertical layers
-      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl), INTENT(in   ) ::   pin    !: input field
-      REAL(wp), DIMENSION(jpi,jpj)           , INTENT(  out) ::   pout   !: output field
-      !
-      INTEGER ::   jk, jl   ! dummy loop indices
-      !!---------------------------------------------------------------------
-      !
-      pout(:,:) = 0._wp
-      DO jl = 1, ksum
-         DO jk = 2, klay 
-            pout(:,:) = pout(:,:) + pin(:,:,jk,jl)
-         END DO
-      END DO
-      !
-   END SUBROUTINE lim_column_sum_energy
-
-
-   SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid )
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE lim_cons_check ***
-      !!
-      !! ** Purpose : Test the conservation of a certain variable
-      !!              For each physical grid cell, check that initial 
-      !!              and final values
-      !!              of a conserved field are equal to within a small value.
-      !!
-      !! ** Method  :
-      !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px1          !: initial field
-      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px2          !: final field
-      REAL(wp)                , INTENT(in   ) ::   pmax_err     !: max allowed error
-      CHARACTER(len=15)       , INTENT(in   ) ::   cd_fieldid   !: field identifyer
-      !
-      INTEGER  ::   ji, jj          ! dummy loop indices
-      INTEGER  ::   inb_error       ! number of g.c where there is a cons. error
-      LOGICAL  ::   llconserv_err   ! = .true. if conservation check failed
-      REAL(wp) ::   zmean_error     ! mean error on error points
-      !!---------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*) ' lim_cons_check '
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
-
-      llconserv_err = .FALSE.
-      inb_error     = 0
-      zmean_error   = 0._wp
-      IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err )   llconserv_err = .TRUE.
-
-      IF( llconserv_err ) THEN
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN
-                  inb_error   = inb_error + 1
-                  zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) )
-                  !
-                  IF(lwp) THEN
-                     WRITE (numout,*) ' ALERTE 99 '
-                     WRITE (numout,*) ' Conservation error: ', cd_fieldid
-                     WRITE (numout,*) ' Point             : ', ji, jj 
-                     WRITE (numout,*) ' lat, lon          : ', gphit(ji,jj), glamt(ji,jj)
-                     WRITE (numout,*) ' Initial value     : ', px1(ji,jj)
-                     WRITE (numout,*) ' Final value       : ', px2(ji,jj)
-                     WRITE (numout,*) ' Difference        : ', px2(ji,jj) - px1(ji,jj)
-                  ENDIF
-               ENDIF
-            END DO
-         END DO
-         !
-      ENDIF
-      IF(lk_mpp)   CALL mpp_sum( inb_error   )
-      IF(lk_mpp)   CALL mpp_sum( zmean_error )
-      !
-      IF( inb_error > 0 .AND. lwp ) THEN
-         zmean_error = zmean_error / REAL( inb_error, wp )
-         WRITE(numout,*) ' Conservation check for : ', cd_fieldid
-         WRITE(numout,*) ' Number of error points : ', inb_error
-         WRITE(numout,*) ' Mean error on these pts: ', zmean_error
-      ENDIF
-      !
-   END SUBROUTINE lim_cons_check
-
 
    SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 8342)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 8355)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
+   SUBROUTINE lim_itd_th_rem( kt )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_th_rem ***
@@ -55,6 +55,4 @@
       !! References : W.H. Lipscomb, JGR 2001
       !!------------------------------------------------------------------
-      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
-      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
       INTEGER , INTENT (in) ::   kt      ! Ocean time step 
       !
@@ -65,5 +63,4 @@
       REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
       REAL(wp) ::   zx3        
-      CHARACTER (len = 15) :: fieldid
 
       INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor   ! donor category index
@@ -75,15 +72,9 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   dummy_es
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume
-      REAL(wp), POINTER, DIMENSION(:)     ::   zvetamin, zvetamax      ! maximum values for etas
       INTEGER , POINTER, DIMENSION(:)     ::   nind_i, nind_j          ! compressed indices for i/j directions
       INTEGER                             ::   nbrem                   ! number of cells with ice to transfer
       REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds"
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zhb0, zhb1              ! category boundaries for thinnes categories
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final   !  ice volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_s_init, vt_s_final   !  snow volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   et_i_init, et_i_final   !  ice energy summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   et_s_init, et_s_final   !  snow energy summed over categories
       INTEGER , POINTER, DIMENSION(:,:)   ::   zremap_flag      ! compute remapping or not ????
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhbnew           ! new boundaries of ice categories
@@ -92,21 +83,9 @@
       CALL wrk_alloc( jpi,jpj, zremap_flag )
       CALL wrk_alloc( jpi,jpj,jpl-1, zdonor )
-      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
+      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b )
       CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )   
       CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
-      CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )   
       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
-      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
-
-      !!----------------------------------------------------------------------------------------------
-      !! 0) Conservation checkand changes in each ice category
-      !!----------------------------------------------------------------------------------------------
-      IF( con_i ) THEN
-         CALL lim_column_sum (jpl,   v_i, vt_i_init)
-         CALL lim_column_sum (jpl,   v_s, vt_s_init)
-         CALL lim_column_sum_energy (jpl, nlay_i,   e_i, et_i_init)
-         dummy_es(:,:,:) = e_s(:,:,1,:)
-         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
-      ENDIF
+      CALL wrk_alloc( jpi,jpj, zhb0, zhb1 )
 
       !!----------------------------------------------------------------------------------------------
@@ -117,10 +96,8 @@
          WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution'
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) ' klbnd :       ', klbnd
-         WRITE(numout,*) ' kubnd :       ', kubnd
       ENDIF
 
       zdhice(:,:,:) = 0._wp
-      DO jl = klbnd, kubnd
+      DO jl = 1, jpl
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -134,11 +111,11 @@
       END DO
 
-      !-----------------------------------------------------------------------------------------------
-      !  2) Compute fractional ice area in each grid cell
-      !-----------------------------------------------------------------------------------------------
-      at_i(:,:) = 0._wp
-      DO jl = klbnd, kubnd
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
+!!      !-----------------------------------------------------------------------------------------------
+!!      !  2) Compute fractional ice area in each grid cell
+!!      !-----------------------------------------------------------------------------------------------
+!!      at_i(:,:) = 0._wp
+!!      DO jl = 1, jpl
+!!         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+!!      END DO
 
       !-----------------------------------------------------------------------------------------------
@@ -160,50 +137,44 @@
 
       !-----------------------------------------------------------------------------------------------
-      !  4) Compute new category boundaries
-      !-----------------------------------------------------------------------------------------------
-      !- 4.1 Compute category boundaries
+      !  4) Compute new category boundaries: 1:jpl-1
+      !-----------------------------------------------------------------------------------------------
       zhbnew(:,:,:) = 0._wp
 
-      DO jl = klbnd, kubnd - 1
-         DO ji = 1, nbrem
-            ii = nind_i(ji)
-            ij = nind_j(ji)
-            !
-            zhbnew(ii,ij,jl) = hi_max(jl)
-            IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN
-               !interpolate between adjacent category growth rates
-               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
-               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
-            ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN
-               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
-            ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN
-               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
-            ENDIF
-         END DO
-
-         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
-         DO ji = 1, nbrem
-            ii = nind_i(ji)
-            ij = nind_j(ji)
-
-            ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible 
-            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
-            IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN
-               zremap_flag(ii,ij) = 0
-            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN
-               zremap_flag(ii,ij) = 0
-            ENDIF
-
-            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max  
-            IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
-            IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0
-            ! clem bug: why is not the following instead?
-            !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
-            !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0
- 
-         END DO
-
-      END DO
-
+      IF( nbrem > 0 ) THEN
+         DO jl = 1, jpl - 1
+            DO ji = 1, nbrem
+               ii = nind_i(ji)
+               ij = nind_j(ji)
+               !
+               ! --- New boundary categories Hn* = Hn + Fn*dt ---
+               !!    Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b)
+               IF    ( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl+1) & a_i_b(jl) /= 0
+                  zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
+                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
+!!               ELSEIF( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,jj,jl+1) <= epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt
+               ELSEIF( a_i_b(ii,ij,jl) >  epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt
+                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
+!!               ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt
+               ELSEIF( a_i_b(ii,ij,jl+1) > epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt
+                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
+               ELSE                                                                       !! a_i_b(jl+1) & a_i_b(jl) = 0 
+                  zhbnew(ii,ij,jl) = hi_max(jl)
+               ENDIF
+            
+               ! --- 2 conditions for remapping --- !
+               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi               
+               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
+               !          in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
+               IF( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) )   zremap_flag(ii,ij) = 0
+               IF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) )   zremap_flag(ii,ij) = 0
+               
+               ! 2) Hn-1 < Hn* < Hn+1  
+               IF( zhbnew(ii,ij,jl) < hi_max(jl-1) )   zremap_flag(ii,ij) = 0
+               IF( zhbnew(ii,ij,jl) > hi_max(jl+1) )   zremap_flag(ii,ij) = 0
+               
+            END DO
+         END DO
+      ENDIF
+   
       !-----------------------------------------------------------------------------------------------
       !  5) Identify cells where ITD is to be remapped
@@ -221,5 +192,5 @@
 
       !-----------------------------------------------------------------------------------------------
-      !  6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
+      !  6) Compute new category boundaries: jpl
       !-----------------------------------------------------------------------------------------------
       DO jj = 1, jpj
@@ -228,18 +199,18 @@
             zhb1(ji,jj) = hi_max(1)
 
-            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN
-               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) )
+            ! boundary for jpl
+            IF( a_i(ji,jj,jpl) > epsi10 ) THEN
+               zhbnew(ji,jj,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ji,jj,jpl) - 2._wp * zhbnew(ji,jj,jpl-1) )
             ELSE
-!clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd)  
-               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway
-            ENDIF
-
-            ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible 
-            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
-            IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN
-               zremap_flag(ji,jj) = 0
-            ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN
-               zremap_flag(ji,jj) = 0
-            ENDIF
+!clem bug               zhbnew(ji,jj,jpl) = hi_max(jpl)  
+               zhbnew(ji,jj,jpl) = hi_max(jpl-1) ! not used anyway
+            ENDIF
+
+            ! 1 condition for remapping for the 1st category
+            ! H0+epsi < h1(t) < H1-epsi 
+            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
+            !    in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
+            IF( zht_i_b(ji,jj,1) < ( zhb0(ji,jj) + epsi10 ) )   zremap_flag(ji,jj) = 0
+            IF( zht_i_b(ji,jj,1) > ( zhb1(ji,jj) - epsi10 ) )   zremap_flag(ji,jj) = 0
 
          END DO
@@ -250,37 +221,37 @@
       !-----------------------------------------------------------------------------------------------
       !- 7.1 g(h) for category 1 at start of time step
-      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   &
-         &                  hR(:,:,klbnd), zremap_flag )
-
-      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd)
+      CALL lim_itd_fitline( 1, zhb0, zhb1, zht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1),   &
+         &                  hR(:,:,1), zremap_flag )
+
+      !- 7.2 Area lost due to melting of thin ice (first category,  1)
       DO ji = 1, nbrem
          ii = nind_i(ji) 
          ij = nind_j(ji) 
 
-         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
-
-            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
+         IF( a_i(ii,ij,1) > epsi10 ) THEN
+
+            zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category
 
             IF( zdh0 < 0.0 ) THEN      !remove area from category 1
-               zdh0 = MIN( -zdh0, hi_max(klbnd) )
+               zdh0 = MIN( -zdh0, hi_max(1) )
                !Integrate g(1) from 0 to dh0 to estimate area melted
-               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd)
+               zetamax = MIN( zdh0, hR(ii,ij,1) ) - hL(ii,ij,1)
 
                IF( zetamax > 0.0 ) THEN
                   zx1    = zetamax
                   zx2    = 0.5 * zetamax * zetamax 
-                  zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed
-                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i                
+                  zda0   = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1                        ! ice area removed
+                  zdamax = a_i(ii,ij,1) * (1.0 - ht_i(ii,ij,1) / zht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i                
                   zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting 
                                                                                                 !     of thin ice (zdamax > 0)
                   ! Remove area, conserving volume
-                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
-                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0
-                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
+                  ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1) - zda0 )
+                  a_i(ii,ij,1)  = a_i(ii,ij,1) - zda0
+                  v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clem-useless ?
                ENDIF
 
             ELSE ! if ice accretion zdh0 > 0
                ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
-               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 
+               zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) ) 
             ENDIF
 
@@ -290,5 +261,5 @@
 
       !- 7.3 g(h) for each thickness category  
-      DO jl = klbnd, kubnd
+      DO jl = 1, jpl
          CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  &
             &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
@@ -296,8 +267,8 @@
 
       !-----------------------------------------------------------------------------------------------
-      !  8) Compute area and volume to be shifted across each boundary
-      !-----------------------------------------------------------------------------------------------
-
-      DO jl = klbnd, kubnd - 1
+      !  8) Compute area and volume to be shifted across each boundary (Eq. 18)
+      !-----------------------------------------------------------------------------------------------
+
+      DO jl = 1, jpl - 1
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -312,20 +283,15 @@
             ij = nind_j(ji)
 
-            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1
-               ! left and right integration limits in eta space
-               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)
-               zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)
+            ! left and right integration limits in eta space
+            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
+               zetamin = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)       ! hi_max(jl) - hL 
+               zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) ! hR - hL
                zdonor(ii,ij,jl) = jl
-
-            ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
-               ! left and right integration limits in eta space
-               zvetamin(ji) = 0.0
-               zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)
+            ELSE                                    ! Hn* <= Hn => transfer from jl+1 to jl
+               zetamin = 0.0
+               zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)   ! hi_max(jl) - hL
                zdonor(ii,ij,jl) = jl + 1
-
-            ENDIF
-
-            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin
-            zetamin = zvetamin(ji)
+            ENDIF
+            zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
 
             zx1  = zetamax - zetamin
@@ -346,5 +312,5 @@
       !! 9) Shift ice between categories
       !!----------------------------------------------------------------------------------------------
-      CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
+      CALL lim_itd_shiftice ( zdonor, zdaice, zdvice )
 
       !!----------------------------------------------------------------------------------------------
@@ -366,34 +332,11 @@
       END DO
 
-      !!----------------------------------------------------------------------------------------------
-      !! 11) Conservation check
-      !!----------------------------------------------------------------------------------------------
-      IF ( con_i ) THEN
-         CALL lim_column_sum (jpl,   v_i, vt_i_final)
-         fieldid = ' v_i : limitd_th '
-         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
-
-         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, et_i_final)
-         fieldid = ' e_i : limitd_th '
-         CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid) 
-
-         CALL lim_column_sum (jpl,   v_s, vt_s_final)
-         fieldid = ' v_s : limitd_th '
-         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
-
-         dummy_es(:,:,:) = e_s(:,:,1,:)
-         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
-         fieldid = ' e_s : limitd_th '
-         CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 
-      ENDIF
-
       CALL wrk_dealloc( jpi,jpj, zremap_flag )
       CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )
-      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
+      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b )
       CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )   
       CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
-      CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )   
       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
-      CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
+      CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 )
 
    END SUBROUTINE lim_itd_th_rem
@@ -404,10 +347,9 @@
       !!                ***  ROUTINE lim_itd_fitline ***
       !!
-      !! ** Purpose :   fit g(h) with a line using area, volume constraints
-      !!
-      !! ** Method  :   Fit g(h) with a line, satisfying area and volume constraints.
-      !!              To reduce roundoff errors caused by large values of g0 and g1,
-      !!              we actually compute g(eta), where eta = h - hL, and hL is the
-      !!              left boundary.
+      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9)
+      !!
+      !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0
+      !!                with eta = h - HL
+      !!                
       !!------------------------------------------------------------------
       INTEGER                     , INTENT(in   ) ::   num_cat      ! category index
@@ -436,10 +378,11 @@
                hR(ji,jj) = HbR(ji,jj)
 
-               ! Change hL or hR if hice falls outside central third of range
-               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )
-               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )
-
-               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj)
-               ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj)
+               ! Change hL or hR if hice falls outside central third of range,
+               ! so that hice is in the central third of the range [HL HR]
+               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) +       hR(ji,jj) )
+               zh23 = 1.0 / 3.0 * (       hL(ji,jj) + 2.0 * hR(ji,jj) )
+
+               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) ! move HR to the left
+               ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) ! move HL to the right
                ENDIF
 
@@ -448,6 +391,6 @@
                zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr
                zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr
-               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )
-               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 )
+               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14
+               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14
                !
             ELSE  ! remap_flag = .false. or a_i < epsi10 
@@ -464,5 +407,5 @@
 
 
-   SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
+   SUBROUTINE lim_itd_shiftice( zdonor, zdaice, zdvice )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_shiftice ***
@@ -473,6 +416,4 @@
       !! ** Method  :
       !!------------------------------------------------------------------
-      INTEGER                           , INTENT(in   ) ::   klbnd    ! Start thickness category index point
-      INTEGER                           , INTENT(in   ) ::   kubnd    ! End point on which the  the computation is applied
       INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in   ) ::   zdonor   ! donor category index
       REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdaice   ! ice area transferred across boundary
@@ -508,5 +449,5 @@
       !----------------------------------------------------------------------------------------------
 
-      DO jl = klbnd, kubnd
+      DO jl = 1, jpl
          zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
       END DO
@@ -516,5 +457,5 @@
       !-------------------------------------------------------------------------------
 
-      DO jl = klbnd, kubnd - 1
+      DO jl = 1, jpl - 1
          nbrem = 0
          DO jj = 1, jpj
@@ -539,47 +480,33 @@
             ENDIF
 
-            !--------------
             ! Ice areas
-            !--------------
             a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl)
             a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl)
 
-            !--------------
             ! Ice volumes
-            !--------------
             v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 
             v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl)
 
-            !--------------
             ! Snow volumes
-            !--------------
             zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij)
             v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow
             v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 
 
-            !--------------------
             ! Snow heat content  
-            !--------------------
             zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij)
             e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow
             e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow
 
-            !--------------
             ! Ice age 
-            !--------------
             zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl)
             oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice
             oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice
 
-            !--------------
             ! Ice salinity
-            !--------------
             zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij)
             smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice
             smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice
 
-            !---------------------
             ! Surface temperature
-            !---------------------
             zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl)
             zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf
@@ -588,14 +515,10 @@
             ! MV MP 2016 
             IF ( nn_pnd_scheme > 0 ) THEN
-            !---------------------
             ! Pond fraction
-            !---------------------
                zdapnd             = a_ip(ii,ij,jl1) * zdaice(ii,ij,jl)
                a_ip(ii,ij,jl1)    = a_ip(ii,ij,jl1) - zdapnd
                a_ip(ii,ij,jl2)    = a_ip(ii,ij,jl2) + zdapnd
 
-            !---------------------
             ! Pond volume
-            !---------------------
                zdvpnd             = v_ip(ii,ij,jl1) * zdaice(ii,ij,jl)
                v_ip(ii,ij,jl1)    = v_ip(ii,ij,jl1) - zdvpnd
@@ -607,8 +530,5 @@
          END DO
 
-         !------------------
          ! Ice heat content
-         !------------------
-
          DO jk = 1, nlay_i
             DO ji = 1, nbrem
@@ -635,5 +555,5 @@
       !-----------------------------------------------------------------
 
-      DO jl = klbnd, kubnd
+      DO jl = 1, jpl
          DO jj = 1, jpj
             DO ji = 1, jpi 
@@ -656,5 +576,5 @@
    
 
-   SUBROUTINE lim_itd_th_reb( klbnd, kubnd )
+   SUBROUTINE lim_itd_th_reb
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_th_reb ***
@@ -664,32 +584,18 @@
       !! ** Method  :
       !!------------------------------------------------------------------
-      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
-      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
-      !
       INTEGER ::   ji,jj, jl   ! dummy loop indices
       INTEGER ::   zshiftflag          ! = .true. if ice must be shifted
-      CHARACTER (len = 15) :: fieldid
 
       INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor           ! donor category index
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice   ! ice area and volume transferred
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories
       !!------------------------------------------------------------------
       
       CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger
       CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice )
-      CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
       !     
-      IF( con_i ) THEN                 ! conservation check
-         CALL lim_column_sum (jpl,   v_i, vt_i_init)
-         CALL lim_column_sum (jpl,   v_s, vt_s_init)
-      ENDIF
-
-      !
       !------------------------------------------------------------------------------
       ! 1) Compute ice thickness.
       !------------------------------------------------------------------------------
-      DO jl = klbnd, kubnd
+      DO jl = 1, jpl
          DO jj = 1, jpj
             DO ji = 1, jpi 
@@ -707,5 +613,5 @@
       ! Initialize shift arrays
       !-------------------------
-      DO jl = klbnd, kubnd
+      DO jl = 1, jpl
          zdonor(:,:,jl) = 0
          zdaice(:,:,jl) = 0._wp
@@ -717,5 +623,5 @@
       !-------------------------
 
-      DO jl = klbnd, kubnd - 1  ! loop over category boundaries
+      DO jl = 1, jpl - 1  ! loop over category boundaries
 
          !---------------------------------------
@@ -742,5 +648,5 @@
 
          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
-            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
+            CALL lim_itd_shiftice( zdonor, zdaice, zdvice )
             ! Reset shift parameters
             zdonor(:,:,jl) = 0
@@ -755,5 +661,5 @@
       !----------------------------
 
-      DO jl = kubnd - 1, 1, -1       ! loop over category boundaries
+      DO jl = jpl - 1, 1, -1       ! loop over category boundaries
 
          !-----------------------------------------
@@ -777,5 +683,5 @@
          
          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
-            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
+            CALL lim_itd_shiftice( zdonor, zdaice, zdvice )
             ! Reset shift parameters
             zdonor(:,:,jl) = 0
@@ -785,22 +691,7 @@
 
       END DO
-
-      !------------------------------------------------------------------------------
-      ! 3) Conservation check
-      !------------------------------------------------------------------------------
-
-      IF( con_i ) THEN
-         CALL lim_column_sum (jpl,   v_i, vt_i_final)
-         fieldid = ' v_i : limitd_reb '
-         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
-
-         CALL lim_column_sum (jpl,   v_s, vt_s_final)
-         fieldid = ' v_s : limitd_reb '
-         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
-      ENDIF
       !
       CALL wrk_dealloc( jpi,jpj,jpl, zdonor )
       CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice )
-      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
 
    END SUBROUTINE lim_itd_th_reb
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8342)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8355)
@@ -319,5 +319,5 @@
       IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt )
+      IF( jpl > 1 )      CALL lim_itd_th_rem( kt )
 
       IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 8342)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 8355)
@@ -104,5 +104,5 @@
       ! Rebin categories with thickness out of bounds
       !----------------------------------------------------
-      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)
+      IF ( jpl > 1 ) CALL lim_itd_th_reb
 
       !-----------------
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 8342)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 8355)
@@ -118,5 +118,5 @@
       ! Rebin categories with thickness out of bounds
       !----------------------------------------------------
-      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
+      IF ( jpl > 1 ) CALL lim_itd_th_reb
 
       !-----------------
