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 8360)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 8361)
@@ -15,5 +15,5 @@
    !!   lim_itd_th_rem   :
    !!   lim_itd_th_reb   :
-   !!   lim_itd_fitline  :
+   !!   lim_itd_glinear  :
    !!   lim_itd_shiftice :
    !!----------------------------------------------------------------------
@@ -24,4 +24,6 @@
    USE ice              ! LIM-3 variables
    USE limvar           ! LIM-3 variables
+   USE limcons          ! conservation tests
+   !
    USE prtctl           ! Print control
    USE in_out_manager   ! I/O manager
@@ -29,5 +31,4 @@
    USE wrk_nemo         ! work arrays
    USE lib_fortran      ! to use key_nosignedzero
-   USE limcons          ! conservation tests
 
    IMPLICIT NONE
@@ -71,8 +72,6 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   hL          ! left boundary for the ITD for each thickness
       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(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume
-      INTEGER , POINTER, DIMENSION(:)     ::   nind_i, nind_j          ! compressed indices for i/j directions
-      INTEGER                             ::   nbrem                   ! number of cells with ice to transfer
+      INTEGER , DIMENSION(jpij)  ::   nind_i, nind_j          ! compressed indices for i/j directions
       REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds"
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zhb0, zhb1              ! category boundaries for thinnes categories
@@ -83,13 +82,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 )
+      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR )
       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), nind_i, nind_j ) 
       CALL wrk_alloc( jpi,jpj, zhb0, zhb1 )
 
-      !!----------------------------------------------------------------------------------------------
-      !! 1) Compute thickness and changes in each ice category
-      !!----------------------------------------------------------------------------------------------
       IF( kt == nit000 .AND. lwp) THEN
          WRITE(numout,*)
@@ -98,35 +93,15 @@
       ENDIF
 
-      zdhice(:,:,:) = 0._wp
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               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           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )
-               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: useless IF statement? 
-            END DO
-         END DO
-      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
-
       !-----------------------------------------------------------------------------------------------
       !  3) Identify grid cells with ice
       !-----------------------------------------------------------------------------------------------
-      nbrem = 0
+      nidx = 0 ; idxice(:) = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
             IF ( at_i(ji,jj) > epsi10 ) THEN
-               nbrem         = nbrem + 1
-               nind_i(nbrem) = ji
-               nind_j(nbrem) = jj
+               nidx         = nidx + 1
+               nind_i(nidx) = ji
+               nind_j(nidx) = jj
+               idxice( nidx ) = (jj - 1) * jpi + ji
                zremap_flag(ji,jj) = 1
             ELSE
@@ -139,24 +114,36 @@
       !  4) Compute new category boundaries: 1:jpl-1
       !-----------------------------------------------------------------------------------------------
+      zdhice(:,:,:) = 0._wp
       zhbnew(:,:,:) = 0._wp
-
-      IF( nbrem > 0 ) THEN
+      zhb0(:,:) = hi_max(0)
+      zhb1(:,:) = hi_max(1)
+
+      IF( nidx > 0 ) THEN
+
+         ! Compute thickness change in each ice category
+         DO jl = 1, jpl
+            DO ji = 1, nidx
+               ii = nind_i(ji)
+               ij = nind_j(ji)
+               zdhice(ii,ij,jl) = ht_i(ii,ij,jl) - ht_i_b(ii,ij,jl)
+            END DO
+         END DO
+         
          DO jl = 1, jpl - 1
-            DO ji = 1, nbrem
+            DO ji = 1, nidx
                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
+               ! --- 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) ) / ( ht_i_b(ii,ij,jl+1) - ht_i_b(ii,ij,jl) )
+                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - ht_i_b(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+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
+               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
                   zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
-               ELSE                                                                       !! a_i_b(jl+1) & a_i_b(jl) = 0 
+               ELSE                                                                         ! a_i_b(jl+1) & a_i_b(jl) = 0 
                   zhbnew(ii,ij,jl) = hi_max(jl)
                ENDIF
@@ -165,5 +152,5 @@
                ! 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)
+               !          in lim_itd_glinear 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
@@ -175,4 +162,23 @@
             END DO
          END DO
+         
+         ! --- New boundary for category jpl --- !
+         DO ji = 1, nidx
+            ii = nind_i(ji)
+            ij = nind_j(ji)
+            IF( a_i(ii,ij,jpl) > epsi10 ) THEN
+               zhbnew(ii,ij,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ii,ij,jpl) - 2._wp * zhbnew(ii,ij,jpl-1) )
+            ELSE
+               zhbnew(ii,ij,jpl) = hi_max(jpl)  
+            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_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
+            IF( ht_i_b(ii,ij,1) < ( hi_max(0) + epsi10 ) )   zremap_flag(ii,ij) = 0
+            IF( ht_i_b(ii,ij,1) > ( hi_max(1) - epsi10 ) )   zremap_flag(ii,ij) = 0
+         END DO
+         
       ENDIF
    
@@ -180,11 +186,11 @@
       !  5) Identify cells where ITD is to be remapped
       !-----------------------------------------------------------------------------------------------
-      nbrem = 0
+      nidx = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
             IF( zremap_flag(ji,jj) == 1 ) THEN
-               nbrem         = nbrem + 1
-               nind_i(nbrem) = ji
-               nind_j(nbrem) = jj
+               nidx         = nidx + 1
+               nind_i(nidx) = ji
+               nind_j(nidx) = jj
             ENDIF
          END DO 
@@ -192,150 +198,127 @@
 
       !-----------------------------------------------------------------------------------------------
-      !  6) Compute new category boundaries: jpl
-      !-----------------------------------------------------------------------------------------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zhb0(ji,jj) = hi_max(0)
-            zhb1(ji,jj) = hi_max(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,jpl) = hi_max(jpl)  
-               zhbnew(ji,jj,jpl) = hi_max(jpl-1) ! not used anyway
+      !  7) Compute g(h) 
+      !-----------------------------------------------------------------------------------------------
+      IF( nidx > 0 ) THEN
+
+         !- 7.1 g(h) for category 1 at start of time step
+         CALL lim_itd_glinear( 1, zhb0, zhb1, ht_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, nidx
+            ii = nind_i(ji) 
+            ij = nind_j(ji) 
+            
+            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(1) )
+                  !Integrate g(1) from 0 to dh0 to estimate area melted
+                  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,1) * zx2 + g0(ii,ij,1) * zx1                        ! ice area removed
+                     zdamax = a_i(ii,ij,1) * (1.0 - ht_i(ii,ij,1) / ht_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,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,0) = MIN( zdh0, hi_max(1) ) 
+               ENDIF
+               
             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
-      END DO
-
-      !-----------------------------------------------------------------------------------------------
-      !  7) Compute g(h) 
-      !-----------------------------------------------------------------------------------------------
-      !- 7.1 g(h) for category 1 at start of time step
-      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,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(1) )
-               !Integrate g(1) from 0 to dh0 to estimate area melted
-               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,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,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,0) = MIN( zdh0, hi_max(1) ) 
-            ENDIF
-
-         ENDIF
-
-      END DO
-
-      !- 7.3 g(h) for each thickness category  
-      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 )
-      END DO
-
-      !-----------------------------------------------------------------------------------------------
-      !  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
-               zdonor(ji,jj,jl) = 0
-               zdaice(ji,jj,jl) = 0.0
-               zdvice(ji,jj,jl) = 0.0
-            END DO
-         END DO
-
-         DO ji = 1, nbrem
+            
+         END DO
+         
+         !- 7.3 g(h) for each thickness category  
+         DO jl = 1, jpl
+            CALL lim_itd_glinear( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  &
+               &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
+         END DO
+         
+         !-----------------------------------------------------------------------------------------------
+         !  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
+                  zdonor(ji,jj,jl) = 0
+                  zdaice(ji,jj,jl) = 0.0
+                  zdvice(ji,jj,jl) = 0.0
+               END DO
+            END DO
+            
+            DO ji = 1, nidx
+               ii = nind_i(ji)
+               ij = nind_j(ji)
+               
+               ! 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                                    ! 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( zetamax, zetamin ) ! no transfer if etamax < etamin
+               
+               zx1  = zetamax - zetamin
+               zwk1 = zetamin * zetamin
+               zwk2 = zetamax * zetamax
+               zx2  = 0.5 * ( zwk2 - zwk1 )
+               zwk1 = zwk1 * zetamin
+               zwk2 = zwk2 * zetamax
+               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
+               nd   = zdonor(ii,ij,jl)
+               zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1
+               zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
+               
+            END DO
+         END DO
+         
+         !!----------------------------------------------------------------------------------------------
+         !! 9) Shift ice between categories
+         !!----------------------------------------------------------------------------------------------
+         CALL lim_itd_shiftice ( zdonor, zdaice, zdvice )
+         
+         !!----------------------------------------------------------------------------------------------
+         !! 10) Make sure ht_i >= minimum ice thickness hi_min
+         !!----------------------------------------------------------------------------------------------
+         
+         DO ji = 1, nidx
             ii = nind_i(ji)
             ij = nind_j(ji)
-
-            ! 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                                    ! 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
+            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 
+               ! MV MP 2016
+               IF ( nn_pnd_scheme > 0 ) THEN
+                  a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin
+               ENDIF
+               ! END MV MP 2016
+               ht_i(ii,ij,1) = rn_himin
             ENDIF
-            zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
-
-            zx1  = zetamax - zetamin
-            zwk1 = zetamin * zetamin
-            zwk2 = zetamax * zetamax
-            zx2  = 0.5 * ( zwk2 - zwk1 )
-            zwk1 = zwk1 * zetamin
-            zwk2 = zwk2 * zetamax
-            zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
-            nd   = zdonor(ii,ij,jl)
-            zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1
-            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
-
-         END DO
-      END DO
-
-      !!----------------------------------------------------------------------------------------------
-      !! 9) Shift ice between categories
-      !!----------------------------------------------------------------------------------------------
-      CALL lim_itd_shiftice ( zdonor, zdaice, zdvice )
-
-      !!----------------------------------------------------------------------------------------------
-      !! 10) Make sure ht_i >= minimum ice thickness hi_min
-      !!----------------------------------------------------------------------------------------------
-
-      DO ji = 1, nbrem
-         ii = nind_i(ji)
-         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 
-            ! MV MP 2016
-            IF ( nn_pnd_scheme > 0 ) THEN
-               a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin
-            ENDIF
-            ! END MV MP 2016
-            ht_i(ii,ij,1) = rn_himin
-         ENDIF
-      END DO
+         END DO
+         
+      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 )
+      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR )
       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), nind_i, nind_j ) 
       CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 )
 
@@ -343,7 +326,7 @@
 
 
-   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_itd_fitline ***
+   SUBROUTINE lim_itd_glinear( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE lim_itd_glinear ***
       !!
       !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9)
@@ -404,5 +387,5 @@
       END DO
       !
-   END SUBROUTINE lim_itd_fitline
+   END SUBROUTINE lim_itd_glinear
 
 
@@ -421,141 +404,110 @@
 
       INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices
-      INTEGER ::   ii, ij                     ! indices when changing from 2D-1D is done
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaTsfn
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka            ! temporary array used here
-
-      REAL(wp) ::   zdvsnow, zdesnow   ! snow volume and energy transferred
-      REAL(wp) ::   zdeice             ! ice energy transferred
-      REAL(wp) ::   zdsm_vice          ! ice salinity times volume transferred
-      REAL(wp) ::   zdo_aice           ! ice age times volume transferred
-      REAL(wp) ::   zdaTsf             ! aicen*Tsfcn transferred
-      ! MV MP 2016 
-      REAL(wp) ::   zdapnd             ! pond fraction transferred
-      REAL(wp) ::   zdvpnd             ! pond volume transferred
-      ! END MV MP 2016 
-
-      INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions
-
-      INTEGER  ::   nbrem             ! number of cells with ice to transfer
-      !!------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi,jpj,jpl, zaTsfn )
-      CALL wrk_alloc( jpi,jpj, zworka )
-      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )
+
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zaTsfn
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zworka            ! temporary array used here
+
+      REAL(wp) ::   ztrans             ! ice/snow transferred
+      !!------------------------------------------------------------------
 
       !----------------------------------------------------------------------------------------------
       ! 1) Define a variable equal to a_i*T_su
       !----------------------------------------------------------------------------------------------
-
       DO jl = 1, jpl
-         zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zaTsfn(ji,jj,jl) = a_i(ji,jj,jl) * t_su(ji,jj,jl)
+            END DO
+         END DO
       END DO
-
+      
       !-------------------------------------------------------------------------------
       ! 2) Transfer volume and energy between categories
       !-------------------------------------------------------------------------------
-
       DO jl = 1, jpl - 1
-         nbrem = 0
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny
-                  nbrem = nbrem + 1
-                  nind_i(nbrem) = ji
-                  nind_j(nbrem) = jj
-               ENDIF
-            END DO
-         END DO
-
-         DO ji = 1, nbrem 
-            ii = nind_i(ji)
-            ij = nind_j(ji)
-
-            jl1 = zdonor(ii,ij,jl)
-            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) )
-            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch
-            IF( jl1 == jl) THEN   ;   jl2 = jl1+1
-            ELSE                  ;   jl2 = jl 
-            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
-            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf 
-
-            ! 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
-               v_ip(ii,ij,jl2)    = v_ip(ii,ij,jl2) + zdvpnd
-
-            ENDIF
-            ! END MV MP 2016 
-
-         END DO
-
+
+               jl1 = zdonor(ji,jj,jl)
+               rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl1) - epsi10 ) )
+               zworka(ji,jj) = zdvice(ji,jj,jl) / MAX( v_i(ji,jj,jl1), epsi10 ) * rswitch
+
+               IF( jl1 == jl) THEN   ;   jl2 = jl1+1
+               ELSE                  ;   jl2 = jl 
+               ENDIF
+               
+               ! Ice areas
+               a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - zdaice(ji,jj,jl)
+               a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + zdaice(ji,jj,jl)
+               
+               ! Ice volumes
+               v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - zdvice(ji,jj,jl) 
+               v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + zdvice(ji,jj,jl)
+               
+               ! Snow volumes
+               ztrans         = v_s(ji,jj,jl1) * zworka(ji,jj)
+               v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - ztrans
+               v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + ztrans 
+               
+               ! Snow heat content  
+               ztrans             = e_s(ji,jj,1,jl1) * zworka(ji,jj)
+               e_s(ji,jj,1,jl1)   = e_s(ji,jj,1,jl1) - ztrans
+               e_s(ji,jj,1,jl2)   = e_s(ji,jj,1,jl2) + ztrans
+               
+               ! Ice age 
+               ztrans             = oa_i(ji,jj,jl1) * zdaice(ji,jj,jl)
+               oa_i(ji,jj,jl1)    = oa_i(ji,jj,jl1) - ztrans
+               oa_i(ji,jj,jl2)    = oa_i(ji,jj,jl2) + ztrans
+               
+               ! Ice salinity
+               ztrans             = smv_i(ji,jj,jl1) * zworka(ji,jj)
+               smv_i(ji,jj,jl1)   = smv_i(ji,jj,jl1) - ztrans
+               smv_i(ji,jj,jl2)   = smv_i(ji,jj,jl2) + ztrans
+               
+               ! Surface temperature
+               ztrans             = t_su(ji,jj,jl1) * zdaice(ji,jj,jl)
+               zaTsfn(ji,jj,jl1)  = zaTsfn(ji,jj,jl1) - ztrans
+               zaTsfn(ji,jj,jl2)  = zaTsfn(ji,jj,jl2) + ztrans
+               
+               ! MV MP 2016 
+               IF ( nn_pnd_scheme > 0 ) THEN
+                  ! Pond fraction
+                  ztrans          = a_ip(ji,jj,jl1) * zdaice(ji,jj,jl)
+                  a_ip(ji,jj,jl1) = a_ip(ji,jj,jl1) - ztrans
+                  a_ip(ji,jj,jl2) = a_ip(ji,jj,jl2) + ztrans
+                  
+                  ! Pond volume
+                  ztrans          = v_ip(ji,jj,jl1) * zdaice(ji,jj,jl)
+                  v_ip(ji,jj,jl1) = v_ip(ji,jj,jl1) - ztrans
+                  v_ip(ji,jj,jl2) = v_ip(ji,jj,jl2) + ztrans
+               ENDIF
+               ! END MV MP 2016 
+               
+            END DO
+         END DO
+         
          ! Ice heat content
          DO jk = 1, nlay_i
-            DO ji = 1, nbrem
-               ii = nind_i(ji)
-               ij = nind_j(ji)
-
-               jl1 = zdonor(ii,ij,jl)
-               IF (jl1 == jl) THEN
-                  jl2 = jl+1
-               ELSE             ! n1 = n+1
-                  jl2 = jl 
-               ENDIF
-
-               zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij)
-               e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice
-               e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice 
-            END DO
-         END DO
-
-      END DO                   ! boundaries, 1 to ncat-1
-
-      !-----------------------------------------------------------------
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  
+                  jl1 = zdonor(ji,jj,jl)
+                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
+                  ELSE                ;  jl2 = jl
+                  ENDIF
+                  
+                  ztrans            = e_i(ji,jj,jk,jl1) * zworka(ji,jj)
+                  e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - ztrans
+                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + ztrans 
+               END DO
+            END DO
+         END DO
+         
+      END DO                   ! boundaries, 1 to jpl-1
+      
       ! Update ice thickness and temperature
-      !-----------------------------------------------------------------
-
       DO jl = 1, jpl
          DO jj = 1, jpj
-            DO ji = 1, jpi 
+            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) 
@@ -567,9 +519,5 @@
             END DO
          END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn )
-      CALL wrk_dealloc( jpi,jpj, zworka )
-      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )
+      END DO      
       !
    END SUBROUTINE lim_itd_shiftice
@@ -584,14 +532,11 @@
       !! ** Method  :
       !!------------------------------------------------------------------
-      INTEGER ::   ji,jj, jl   ! dummy loop indices
+      INTEGER ::   ji, jj, jl   ! dummy loop indices
       INTEGER ::   zshiftflag          ! = .true. if ice must be shifted
 
-      INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor           ! donor category index
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice   ! ice area and volume transferred
+      INTEGER , DIMENSION(jpi,jpj,jpl) ::   zdonor           ! donor category index
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zdaice, zdvice   ! ice area and volume transferred
       !!------------------------------------------------------------------
       
-      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger
-      CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice )
-      !     
       !------------------------------------------------------------------------------
       ! 1) Compute ice thickness.
@@ -635,8 +580,4 @@
                   zshiftflag        = 1
                   zdonor(ji,jj,jl)  = jl 
-                  ! begin TECLIM change
-                  !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) + epsi20 ) / ht_i(ji,jj,jl)  
@@ -645,5 +586,5 @@
             END DO
          END DO
-         IF(lk_mpp)   CALL mpp_max( zshiftflag )
+         IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ???
 
          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
@@ -680,5 +621,5 @@
          END DO
 
-         IF(lk_mpp)   CALL mpp_max( zshiftflag )
+         IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ???
          
          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
@@ -692,7 +633,4 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpl, zdonor )
-      CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice )
-
    END SUBROUTINE lim_itd_th_reb
 
@@ -704,6 +642,6 @@
    SUBROUTINE lim_itd_th_rem
    END SUBROUTINE lim_itd_th_rem
-   SUBROUTINE lim_itd_fitline
-   END SUBROUTINE lim_itd_fitline
+   SUBROUTINE lim_itd_glinear
+   END SUBROUTINE lim_itd_glinear
    SUBROUTINE lim_itd_shiftice
    END SUBROUTINE lim_itd_shiftice
