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 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 8369)
@@ -25,4 +25,5 @@
    USE limvar           ! LIM-3 variables
    USE limcons          ! conservation tests
+   USE limtab
    !
    USE prtctl           ! Print control
@@ -35,6 +36,6 @@
    PRIVATE
 
-   PUBLIC   lim_itd_th_rem
-   PUBLIC   lim_itd_th_reb
+   PUBLIC   lim_itd_th_rem   ! called in limthd
+   PUBLIC   lim_itd_th_reb   ! called in limupdate
 
    !!----------------------------------------------------------------------
@@ -58,32 +59,20 @@
       INTEGER , INTENT (in) ::   kt      ! Ocean time step 
       !
-      INTEGER  ::   ji, jj, jl     ! dummy loop index
-      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
-      INTEGER  ::   nd             ! local integer
+      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index
+      INTEGER  ::   nidx2                ! local integer
       REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
       REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
       REAL(wp) ::   zx3        
-
-      INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor   ! donor category index
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdhice      ! ice thickness increment
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   g0          ! coefficients for fitting the line of the ITD
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   g1          ! coefficients for fitting the line of the ITD
-      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(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume
-      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
-      INTEGER , POINTER, DIMENSION(:,:)   ::   zremap_flag      ! compute remapping or not ????
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhbnew           ! new boundaries of ice categories
-      !!------------------------------------------------------------------
-
-      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 )
-      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )   
-      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
-      CALL wrk_alloc( jpi,jpj, zhb0, zhb1 )
+      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
+
+      INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not
+      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
+      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
+      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD
+      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness
+      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume
+      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories
+      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories
+      !!------------------------------------------------------------------
 
       IF( kt == nit000 .AND. lwp) THEN
@@ -94,5 +83,5 @@
 
       !-----------------------------------------------------------------------------------------------
-      !  3) Identify grid cells with ice
+      !  1) Identify grid cells with ice
       !-----------------------------------------------------------------------------------------------
       nidx = 0 ; idxice(:) = 0
@@ -101,50 +90,46 @@
             IF ( at_i(ji,jj) > epsi10 ) THEN
                nidx         = nidx + 1
-               nind_i(nidx) = ji
-               nind_j(nidx) = jj
                idxice( nidx ) = (jj - 1) * jpi + ji
-               zremap_flag(ji,jj) = 1
-            ELSE
-               zremap_flag(ji,jj) = 0
             ENDIF
          END DO
       END DO
-
+      
       !-----------------------------------------------------------------------------------------------
-      !  4) Compute new category boundaries: 1:jpl-1
+      !  2) Compute new category boundaries
       !-----------------------------------------------------------------------------------------------
-      zdhice(:,:,:) = 0._wp
-      zhbnew(:,:,:) = 0._wp
-      zhb0(:,:) = hi_max(0)
-      zhb1(:,:) = hi_max(1)
-
       IF( nidx > 0 ) THEN
 
-         ! Compute thickness change in each ice category
+         zdhice(:,:) = 0._wp
+         zhbnew(:,:) = 0._wp
+
+         CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), ht_ib_2d(1:nidx,1:jpl), ht_i_b )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b  )
+
          DO jl = 1, jpl
+            ! Compute thickness change in each ice category
             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)
+               zdhice(ji,jl) = ht_i_2d(ji,jl) - ht_ib_2d(ji,jl)
             END DO
          END DO
          
+         ! --- New boundaries for category 1:jpl-1 --- !
          DO jl = 1, jpl - 1
+
             DO ji = 1, nidx
-               ii = nind_i(ji)
-               ij = nind_j(ji)
                !
-               ! --- New boundary categories Hn* = Hn + Fn*dt --- !
+               ! --- New boundary: 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
-                  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)
+               IF    ( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl+1) & a(jl) /= 0
+                  zslope           = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( ht_ib_2d(ji,jl+1) - ht_ib_2d(ji,jl) )
+                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - ht_ib_2d(ji,jl) )
+               ELSEIF( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN   ! a(jl+1)=0 => Hn* = Hn + fn*dt
+                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl)
+               ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl)=0 => Hn* = Hn + fn+1*dt
+                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1)
+               ELSE                                                                       ! a(jl+1) & a(jl) = 0 
+                  zhbnew(ji,jl) = hi_max(jl)
                ENDIF
             
@@ -153,125 +138,130 @@
                !    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_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
+               IF( a_i_2d(ji,jl  ) > epsi10 .AND. ht_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   idxice(ji) = 0
+               IF( a_i_2d(ji,jl+1) > epsi10 .AND. ht_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   idxice(ji) = 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
+               IF( zhbnew(ji,jl) < hi_max(jl-1) )   idxice(ji) = 0
+               IF( zhbnew(ji,jl) > hi_max(jl+1) )   idxice(ji) = 0
                
             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) )
+         ! --- New boundaries for category jpl --- !
+         DO ji = 1, nidx
+            IF( a_i_2d(ji,jpl) > epsi10 ) THEN
+               zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
             ELSE
-               zhbnew(ii,ij,jpl) = hi_max(jpl)  
-            ENDIF
-            
-            ! 1 condition for remapping for the 1st category
+               zhbnew(ji,jpl) = hi_max(jpl)  
+            ENDIF
+            
+            ! --- 1 additional condition for remapping (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
+            IF( ht_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   idxice(ji) = 0
+            IF( ht_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   idxice(ji) = 0
+         END DO
+
+         !-----------------------------------------------------------------------------------------------
+         !  3) Identify cells where remapping
+         !-----------------------------------------------------------------------------------------------
+         nidx2 = 0 ; idxice2(:) = 0
+         DO ji = 1, nidx
+            IF( idxice(ji) /= 0 ) THEN
+               nidx2 = nidx2 + 1
+               idxice2(nidx2) = idxice(ji)
+               zhbnew(nidx2,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices
+            ENDIF
+         END DO
+         idxice(:) = idxice2(:)
+         nidx      = nidx2
          
       ENDIF
    
+
       !-----------------------------------------------------------------------------------------------
-      !  5) Identify cells where ITD is to be remapped
-      !-----------------------------------------------------------------------------------------------
-      nidx = 0
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( zremap_flag(ji,jj) == 1 ) THEN
-               nidx         = nidx + 1
-               nind_i(nidx) = ji
-               nind_j(nidx) = jj
-            ENDIF
-         END DO 
-      END DO 
-
-      !-----------------------------------------------------------------------------------------------
-      !  7) Compute g(h) 
+      !  4) 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)
+         zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1)
+         g0(:,:) = 0._wp     ; g1(:,:) = 0._wp 
+         hL(:,:) = 0._wp     ; hR(:,:) = 0._wp 
+         
+         DO jl = 1, jpl
+
+            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) )
+            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
+            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
+            CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
+
+            IF( jl == 1 ) THEN
+               
+               ! --- g(h) for category 1 --- !
+               CALL lim_itd_glinear( zhb0(1:nidx), zhb1(1:nidx), ht_ib_1d(1:nidx), a_i_1d(1:nidx),  &   ! in
+                  &                  g0(1:nidx,1), g1(1:nidx,1), hL(1:nidx,1)    , hR(1:nidx,1) )       ! out
+               
+               ! Area lost due to melting of thin ice
+               DO ji = 1, nidx
                   
-                  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 ?
+                  IF( a_i_1d(ji) > epsi10 ) THEN
+                     
+                     zdh0 =  ht_i_1d(ji) - ht_ib_1d(ji)                
+                     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(ji,1) ) - hL(ji,1)
+                        
+                        IF( zetamax > 0.0 ) THEN
+                           zx1    = zetamax
+                           zx2    = 0.5 * zetamax * zetamax 
+                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
+                           zdamax = a_i_1d(ji) * (1.0 - ht_i_1d(ji) / ht_ib_1d(ji) ) ! 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_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
+                           a_i_1d(ji)  = a_i_1d(ji) - zda0
+                           v_i_1d(ji)  = a_i_1d(ji) * ht_i_1d(ji) ! 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(ji,0) = MIN( zdh0, hi_max(1) ) 
+                     ENDIF
+                     
                   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_glinear( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  &
-               &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
+               END DO
+
+               CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
+               CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
+               CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
+               
+            ENDIF ! jl=1
+            
+            ! --- g(h) for each thickness category --- !  
+            CALL lim_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx), a_i_1d(1:nidx),  &  ! in
+               &                  g0(1:nidx,jl)      , g1(1:nidx,jl)    , hL(1:nidx,jl)  , hR(1:nidx,jl) )     ! out
+
          END DO
          
          !-----------------------------------------------------------------------------------------------
-         !  8) Compute area and volume to be shifted across each boundary (Eq. 18)
+         !  5) 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
+               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
+                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL 
+                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL
+                  jdonor(ji,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
+                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL
+                  jdonor(ji,jl) = jl + 1
                ENDIF
                zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
@@ -284,47 +274,45 @@
                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)
+               jcat = jdonor(ji,jl)
+               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1
+               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat)
                
             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)
-            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 
+         !----------------------------------------------------------------------------------------------
+         ! 6) Shift ice between categories
+         !----------------------------------------------------------------------------------------------
+         CALL lim_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )
+         
+         !----------------------------------------------------------------------------------------------
+         ! 7) Make sure ht_i >= minimum ice thickness hi_min
+         !----------------------------------------------------------------------------------------------
+         CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
+         
+         DO ji = 1, nidx
+            IF ( a_i_1d(ji) > epsi10 .AND. ht_i_1d(ji) < rn_himin ) THEN
+               a_i_1d (ji) = a_i_1d(ji) * ht_i_1d(ji) / 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
+                  a_ip_1d(ji) = a_ip_1d(ji) * ht_i_1d(ji) / rn_himin
                ENDIF
                ! END MV MP 2016
-               ht_i(ii,ij,1) = rn_himin
-            ENDIF
-         END DO
-         
+               ht_i_1d(ji) = rn_himin
+            ENDIF
+         END DO
+
+         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
+        
       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 )
-      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )   
-      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
-      CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 )
-
    END SUBROUTINE lim_itd_th_rem
 
 
-   SUBROUTINE lim_itd_glinear( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
+   SUBROUTINE lim_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_glinear ***
@@ -336,13 +324,10 @@
       !!                
       !!------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   num_cat      ! category index
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   HbL, HbR     ! left and right category boundaries
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   hice         ! ice thickness
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   g0, g1       ! coefficients in linear equation for g(eta)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hL           ! min value of range over which g(h) > 0
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hR           ! max value of range over which g(h) > 0
-      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  !
+      REAL(wp), DIMENSION(:), INTENT(in   ) ::   HbL, HbR      ! left and right category boundaries
+      REAL(wp), DIMENSION(:), INTENT(in   ) ::   phice, paice  ! ice thickness and concentration
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   pg0, pg1      ! coefficients in linear equation for g(eta)
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   phL, phR      ! min and max value of range over which g(h) > 0
       !
-      INTEGER  ::   ji,jj        ! horizontal indices
+      INTEGER  ::   ji           ! horizontal indices
       REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL)
       REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL)
@@ -351,44 +336,41 @@
       !!------------------------------------------------------------------
       !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
+         DO ji = 1, nidx
             !
-            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   &
-               &                        .AND. hice(ji,jj)        > 0._wp )  THEN
+            IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN
 
                ! Initialize hL and hR
-               hL(ji,jj) = HbL(ji,jj)
-               hR(ji,jj) = HbR(ji,jj)
+               phL(ji) = HbL(ji)
+               phR(ji) = HbR(ji)
 
                ! 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
+               zh13 = 1.0 / 3.0 * ( 2.0 * phL(ji) +       phR(ji) )
+               zh23 = 1.0 / 3.0 * (       phL(ji) + 2.0 * phR(ji) )
+
+               IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left
+               ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right
                ENDIF
 
                ! Compute coefficients of g(eta) = g0 + g1*eta
-               zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj))
-               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 )      ! Eq. 14
-               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14
+               zdhr = 1._wp / (phR(ji) - phL(ji))
+               zwk1 = 6._wp * paice(ji) * zdhr
+               zwk2 = ( phice(ji) - phL(ji) ) * zdhr
+               pg0(ji) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14
+               pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14
                !
             ELSE  ! remap_flag = .false. or a_i < epsi10 
-               hL(ji,jj) = 0._wp
-               hR(ji,jj) = 0._wp
-               g0(ji,jj) = 0._wp
-               g1(ji,jj) = 0._wp
+               phL(ji) = 0._wp
+               phR(ji) = 0._wp
+               pg0(ji) = 0._wp
+               pg1(ji) = 0._wp
             ENDIF
             !
          END DO
-      END DO
       !
    END SUBROUTINE lim_itd_glinear
 
 
-   SUBROUTINE lim_itd_shiftice( zdonor, zdaice, zdvice )
+   SUBROUTINE lim_itd_shiftice( kdonor, pdaice, pdvice )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_shiftice ***
@@ -399,15 +381,26 @@
       !! ** Method  :
       !!------------------------------------------------------------------
-      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
-      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdvice   ! ice volume transferred across boundary
-
-      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices
-
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zaTsfn
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zworka            ! temporary array used here
+      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
+
+      INTEGER ::   ii,ij, ji, jj, jl, jl2, jl1, jk   ! dummy loop indices
+
+      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn
+      REAL(wp), DIMENSION(jpij)     ::   zworka            ! temporary array used here
+      REAL(wp), DIMENSION(jpij)     ::   zworkv            ! temporary array used here
 
       REAL(wp) ::   ztrans             ! ice/snow transferred
       !!------------------------------------------------------------------
+         
+      CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   )
+      CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su   )
 
       !----------------------------------------------------------------------------------------------
@@ -415,8 +408,6 @@
       !----------------------------------------------------------------------------------------------
       DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zaTsfn(ji,jj,jl) = a_i(ji,jj,jl) * t_su(ji,jj,jl)
-            END DO
+         DO ji = 1, nidx
+            zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
          END DO
       END DO
@@ -426,79 +417,94 @@
       !-------------------------------------------------------------------------------
       DO jl = 1, jpl - 1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               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 
-               
+         DO ji = 1, nidx
+            
+            jl1 = kdonor(ji,jl)
+            IF    ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
+            ELSE                        ;   jl2 = jl 
+            ENDIF
+
+            rswitch    = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) )
+            zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch
+
+            rswitch    = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) )
+            zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch        
+               
+            ! Ice areas
+            a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)
+            a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
+               
+            ! Ice volumes
+            v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl) 
+            v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
+            
+            ! Snow volumes
+            ztrans         = v_s_2d(ji,jl1) * zworkv(ji)
+            v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
+            v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
+                       
+            ! Ice age 
+            ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
+            oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans
+            oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans
+            
+            ! Ice salinity
+            ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)
+            smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans
+            smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans
+            
+            ! Surface temperature
+            ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
+            zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
+            zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
+               
+            ! MV MP 2016 
+            IF ( nn_pnd_scheme > 0 ) THEN
+               ! Pond fraction
+               ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
+               a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
+               a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
+               
+               ! Pond volume (also proportional to da/a)
+               ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
+               v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
+               v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
+            ENDIF
+            ! END MV MP 2016 
+               
+         END DO
+
+         ! Snow heat content
+         DO jk = 1, nlay_s
+
+            DO ji = 1, nidx
+               ii = MOD( idxice(ji) - 1, jpi ) + 1
+               ij = ( idxice(ji) - 1 ) / jpi + 1
+
+               jl1 = kdonor(ji,jl)
+               IF(jl1 == jl) THEN  ;  jl2 = jl+1
+               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 
             END DO
          END DO
+
          
          ! Ice heat content
          DO jk = 1, nlay_i
-            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
+            DO ji = 1, nidx
+               ii = MOD( idxice(ji) - 1, jpi ) + 1
+               ij = ( idxice(ji) - 1 ) / jpi + 1
+               
+               jl1 = kdonor(ji,jl)
+               IF(jl1 == jl) THEN  ;  jl2 = jl+1
+               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 
             END DO
          END DO
@@ -506,18 +512,29 @@
       END DO                   ! boundaries, 1 to jpl-1
       
-      ! Update ice thickness and temperature
+      !-------------------------------------------------------------------------------
+      ! 3) Update ice thickness and temperature
+      !-------------------------------------------------------------------------------
       DO jl = 1, jpl
-         DO jj = 1, jpj
-            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) 
-                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 
-               ELSE
-                  ht_i(ji,jj,jl)  = 0._wp
-                  t_su(ji,jj,jl)  = rt0
-               ENDIF
-            END DO
-         END DO
-      END DO      
+         DO ji = 1, nidx
+            IF ( a_i_2d(ji,jl) > epsi10 ) THEN 
+               ht_i_2d(ji,jl)  =  v_i_2d   (ji,jl) / a_i_2d(ji,jl) 
+               t_su_2d(ji,jl)  =  zaTsfn(ji,jl) / a_i_2d(ji,jl) 
+            ELSE
+               ht_i_2d(ji,jl)  = 0._wp
+               t_su_2d(ji,jl)  = rt0
+            ENDIF
+         END DO
+      END DO
+      
+      CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   )
+      CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su   )
+
       !
    END SUBROUTINE lim_itd_shiftice
@@ -530,103 +547,93 @@
       !! ** Purpose : rebin - rebins thicknesses into defined categories
       !!
-      !! ** Method  :
+      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
+      !!              or entire (for top to down) area, volume, and energy
+      !!              to the neighboring category
       !!------------------------------------------------------------------
       INTEGER ::   ji, jj, jl   ! dummy loop indices
-      INTEGER ::   zshiftflag          ! = .true. if ice must be shifted
-
-      INTEGER , DIMENSION(jpi,jpj,jpl) ::   zdonor           ! donor category index
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zdaice, zdvice   ! ice area and volume transferred
+      !
+      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
+      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
       !!------------------------------------------------------------------
       
-      !------------------------------------------------------------------------------
-      ! 1) Compute ice thickness.
-      !------------------------------------------------------------------------------
-      DO jl = 1, jpl
+      DO jl = 1, jpl-1
+         jdonor(:,jl) = 0
+         zdaice(:,jl) = 0._wp
+         zdvice(:,jl) = 0._wp
+      END DO
+
+      !---------------------------------------
+      ! identify thicknesses that are too big
+      !---------------------------------------
+      DO jl = 1, jpl - 1
+
+         nidx = 0 ; idxice(:) = 0
          DO jj = 1, jpj
-            DO ji = 1, jpi 
-               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
-               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
-            END DO
-         END DO
-      END DO
-
-      !------------------------------------------------------------------------------
-      ! 2) If a category thickness is not in bounds, shift the
-      ! entire area, volume, and energy to the neighboring category
-      !------------------------------------------------------------------------------
-      !-------------------------
-      ! Initialize shift arrays
-      !-------------------------
-      DO jl = 1, jpl
-         zdonor(:,:,jl) = 0
-         zdaice(:,:,jl) = 0._wp
-         zdvice(:,:,jl) = 0._wp
-      END DO
-
-      !-------------------------
-      ! Move thin categories up
-      !-------------------------
-
-      DO jl = 1, jpl - 1  ! loop over category boundaries
-
-         !---------------------------------------
-         ! identify thicknesses that are too big
-         !---------------------------------------
-         zshiftflag = 0
-
-         DO jj = 1, jpj 
-            DO ji = 1, jpi 
-               IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN 
-                  zshiftflag        = 1
-                  zdonor(ji,jj,jl)  = jl 
-                  ! 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)  
-                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 )
-               ENDIF
-            END DO
-         END DO
-         IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ???
-
-         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
-            CALL lim_itd_shiftice( zdonor, zdaice, zdvice )
+            DO ji = 1, jpi
+               IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
+                  nidx = nidx + 1
+                  idxice( nidx ) = (jj - 1) * jpi + ji                  
+               ENDIF
+            ENDDO
+         ENDDO
+         
+!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
+
+         DO ji = 1, nidx
+            jdonor(ji,jl)  = jl 
+            ! how much of a_i you send in cat sup is somewhat arbitrary
+!!clem: these do not work properly after a restart (I do not know why)
+!!          zdaice(ji,jl)  = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji)  
+!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
+!!clem: these do not work properly after a restart (I do not know why)
+!!            zdaice(ji,jl)  = a_i_1d(ji)
+!!            zdvice(ji,jl)  = v_i_1d(ji)
+!!clem: these are from UCL and work ok
+            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
+            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
+            
+         END DO
+
+         IF( nidx > 0 ) THEN
+            CALL lim_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1
             ! Reset shift parameters
-            zdonor(:,:,jl) = 0
-            zdaice(:,:,jl) = 0._wp
-            zdvice(:,:,jl) = 0._wp
+            jdonor(1:nidx,jl) = 0
+            zdaice(1:nidx,jl) = 0._wp
+            zdvice(1:nidx,jl) = 0._wp
          ENDIF
          !
       END DO
 
-      !----------------------------
-      ! Move thick categories down
-      !----------------------------
-
-      DO jl = jpl - 1, 1, -1       ! loop over category boundaries
-
-         !-----------------------------------------
-         ! Identify thicknesses that are too small
-         !-----------------------------------------
-         zshiftflag = 0
-
+      !-----------------------------------------
+      ! Identify thicknesses that are too small
+      !-----------------------------------------
+      DO jl = jpl - 1, 1, -1
+
+         nidx = 0 ; idxice(:) = 0
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
-                  !
-                  zshiftflag = 1
-                  zdonor(ji,jj,jl) = jl + 1
-                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 
-                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
-               ENDIF
-            END DO
-         END DO
-
-         IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ???
-         
-         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
-            CALL lim_itd_shiftice( zdonor, zdaice, zdvice )
+               IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
+                  nidx = nidx + 1
+                  idxice( nidx ) = (jj - 1) * jpi + ji                  
+               ENDIF
+            ENDDO
+         ENDDO
+
+         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok
+         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok
+         DO ji = 1, nidx
+            jdonor(ji,jl) = jl + 1
+            zdaice(ji,jl) = a_i_1d(ji) 
+            zdvice(ji,jl) = v_i_1d(ji)
+         END DO
+         
+         IF( nidx > 0 ) THEN
+            CALL lim_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl
             ! Reset shift parameters
-            zdonor(:,:,jl) = 0
-            zdaice(:,:,jl) = 0._wp
-            zdvice(:,:,jl) = 0._wp
+            jdonor(1:nidx,jl) = 0
+            zdaice(1:nidx,jl) = 0._wp
+            zdvice(1:nidx,jl) = 0._wp
          ENDIF
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8369)
@@ -110,8 +110,8 @@
 
       IF ( .NOT. ln_pnd ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' Melt ponds are not activated '
-         WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. '
-         WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero '
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' Melt ponds are not activated '
+         IF(lwp) WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. '
+         IF(lwp) WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero '
          ln_pnd_rad    = .FALSE.
          ln_pnd_fw     = .FALSE.
@@ -130,5 +130,5 @@
 
       IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN
-         WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false ' 
+         IF(lwp) WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false ' 
          ln_pnd_fw     = .FALSE.
          IF(lwp) THEN                     ! control print
@@ -138,6 +138,6 @@
 
       IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN
-         WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '
-         WRITE(numout,*) ' Run aborted '
+         IF(lwp) WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '
+         IF(lwp) WRITE(numout,*) ' Run aborted '
          CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' )
       ENDIF
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90	(revision 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90	(revision 8369)
@@ -13,10 +13,13 @@
    USE par_kind
    USE par_oce
+   USE ice, ONLY : jpl
    
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tab_2d_1d   ! called by limthd
-   PUBLIC   tab_1d_2d   ! called by limthd
+   PUBLIC   tab_3d_2d
+   PUBLIC   tab_2d_1d
+   PUBLIC   tab_2d_3d
+   PUBLIC   tab_1d_2d
 
    !!----------------------------------------------------------------------
@@ -26,4 +29,25 @@
    !!----------------------------------------------------------------------
 CONTAINS
+
+
+   SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab1d, tab2d )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tab_2d_1d  ***
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   ndim1d   ! 1d size
+      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind  ! input index
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab2d    ! input 2D field
+      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(  out) ::   tab1d    ! output 1D field
+      !
+      INTEGER ::   jl, jn, jid, jjd
+      !!----------------------------------------------------------------------
+      DO jl = 1, jpl
+         DO jn = 1, ndim1d
+            jid          = MOD( tab_ind(jn) - 1 , jpi ) + 1
+            jjd          =    ( tab_ind(jn) - 1 ) / jpi + 1
+            tab1d(jn,jl) = tab2d(jid,jjd,jl)
+         END DO
+      END DO
+   END SUBROUTINE tab_3d_2d
 
    SUBROUTINE tab_2d_1d( ndim1d, tab_ind, tab1d, tab2d )
@@ -45,4 +69,23 @@
    END SUBROUTINE tab_2d_1d
 
+   SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab1d, tab2d )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tab_2d_1d  ***
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   ndim1d    ! 1D size
+      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind   ! input index
+      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab1d     ! input 1D field
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   tab2d     ! output 2D field
+      !
+      INTEGER ::   jl, jn, jid, jjd
+      !!----------------------------------------------------------------------
+      DO jl = 1, jpl
+         DO jn = 1, ndim1d
+            jid               = MOD( tab_ind(jn) - 1 ,  jpi ) + 1
+            jjd               =    ( tab_ind(jn) - 1 ) / jpi  + 1
+            tab2d(jid,jjd,jl) = tab1d(jn,jl)
+         END DO
+      END DO
+   END SUBROUTINE tab_2d_3d
 
    SUBROUTINE tab_1d_2d( ndim1d, tab_ind, tab1d, tab2d )
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 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8369)
@@ -237,41 +237,41 @@
          IF( nidx > 0 ) THEN  ! If there is no ice, do nothing.
             !                                                                
-            s_i_new   (:) = 0._wp ; dh_s_tot (:) = 0._wp                     ! --- some init --- !
+            s_i_new   (:) = 0._wp ; dh_s_tot (:) = 0._wp            ! --- some init --- !
             dh_i_surf (:) = 0._wp ; dh_i_bott(:) = 0._wp
             dh_snowice(:) = 0._wp ; dh_i_sub (:) = 0._wp
 
-                              CALL lim_thd_1d2d( jl, 1 )               ! --- Move to 1D arrays --- !
-            !
-            DO jk = 1, nlay_i                                                ! --- Change units from J/m2 to J/m3 --- !
-               WHERE( ht_i_1d(:)>0._wp ) e_i_1d(:,jk) = e_i_1d(:,jk) / (ht_i_1d(:) * a_i_1d(:)) * nlay_i
+                              CALL lim_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- !
+            !
+            DO jk = 1, nlay_i                                       ! --- Change units from J/m2 to J/m3 --- !
+               WHERE( ht_i_1d(1:nidx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (ht_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i
             ENDDO
             DO jk = 1, nlay_s
-               WHERE( ht_s_1d(:)>0._wp ) e_s_1d(:,jk) = e_s_1d(:,jk) / (ht_s_1d(:) * a_i_1d(:)) * nlay_s
+               WHERE( ht_s_1d(1:nidx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (ht_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s
             ENDDO
             !
-            IF( ln_limdH )    CALL lim_thd_dif                    ! --- Ice/Snow Temperature profile --- !
-            !
-            IF( ln_limdH )    CALL lim_thd_dh                     ! --- Ice/Snow thickness --- !    
+            IF( ln_limdH )    CALL lim_thd_dif                      ! --- Ice/Snow Temperature profile --- !
+            !
+            IF( ln_limdH )    CALL lim_thd_dh                       ! --- Ice/Snow thickness --- !    
             !
             IF( ln_limdH )    CALL lim_thd_ent( e_i_1d(1:nidx,:) )  ! --- Ice enthalpy remapping --- !
             !
-                              CALL lim_thd_sal                    ! --- Ice salinity --- !    
-            !
-                              CALL lim_thd_temp                   ! --- temperature update --- !
+                              CALL lim_thd_sal                      ! --- Ice salinity --- !    
+            !
+                              CALL lim_thd_temp                     ! --- temperature update --- !
             !
             IF( ln_limdH ) THEN
                IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
-                              CALL lim_thd_lam                    ! --- extra lateral melting if monocat --- !
+                              CALL lim_thd_lam                      ! --- extra lateral melting if monocat --- !
                END IF
             END IF
             !
-            DO jk = 1, nlay_i                                                ! --- Change units from J/m3 to J/m2 --- !
-               e_i_1d(:,jk) = e_i_1d(:,jk) * ht_i_1d(:) * a_i_1d(:) * r1_nlay_i
+            DO jk = 1, nlay_i                                       ! --- Change units from J/m3 to J/m2 --- !
+               e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) * ht_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i
             ENDDO
             DO jk = 1, nlay_s
-               e_s_1d(:,jk) = e_s_1d(:,jk) * ht_s_1d(:) * a_i_1d(:) * r1_nlay_s
+               e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) * ht_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s
             ENDDO
             !
-                              CALL lim_thd_1d2d( jl, 2 )               ! --- Move to 2D arrays --- !
+                              CALL lim_thd_1d2d( jl, 2 )            ! --- Move to 2D arrays --- !
             !
             IF( lk_mpp )      CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
@@ -288,5 +288,5 @@
 
       IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_thd_da', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      IF( ln_limdA)           CALL lim_thd_da                                ! --- lateral melting --- !
+      IF( ln_limdA)           CALL lim_thd_da                       ! --- lateral melting --- !
 
       ! Change thickness to volume
@@ -441,5 +441,4 @@
          !
          CALL tab_2d_1d( nidx, idxice(1:nidx), qprec_ice_1d(1:nidx), qprec_ice        )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qevap_ice_1d(1:nidx), qevap_ice(:,:,jl)  )
          CALL tab_2d_1d( nidx, idxice(1:nidx), qsr_ice_1d  (1:nidx), qsr_ice(:,:,jl)  )
          CALL tab_2d_1d( nidx, idxice(1:nidx), fr1_i0_1d   (1:nidx), fr1_i0           )
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90	(revision 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90	(revision 8369)
@@ -163,31 +163,33 @@
             
             DO ji = 1, nidx
-               ! decrease of concentration for the category jl
-               !    each category contributes to melting in proportion to its concentration
-               zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji)
-               
-               ! Contribution to salt flux
-               sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice
-               
-               ! Contribution to heat flux into the ocean [W.m-2], (<0)  
-               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice 
-              
-               ! Contribution to mass flux
-               wfx_lam_1d(ji) =  wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) )
-           
-               !! adjust e_i ???
-               e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
-               e_s_1d(ji,1)        = e_s_1d(ji,1)        * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
- 
-               ! new concentration
-               a_i_1d(ji) = a_i_1d(ji) + zda
-
-               ! ensure that ht_i = 0 where a_i = 0
-               IF( a_i_1d(ji) == 0._wp ) THEN
-                  ht_i_1d(ji) = 0._wp
-                  ht_s_1d(ji) = 0._wp
-               ENDIF            
-            END DO
-
+               IF( a_i_1d(ji) > 0._wp ) THEN
+                  ! decrease of concentration for the category jl
+                  !    each category contributes to melting in proportion to its concentration
+                  zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji)
+                  
+                  ! Contribution to salt flux
+                  sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice
+                  
+                  ! Contribution to heat flux into the ocean [W.m-2], (<0)  
+                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice 
+                  
+                  ! Contribution to mass flux
+                  wfx_lam_1d(ji) =  wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) )
+                  
+                  !! adjust e_i ???
+                  e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
+                  e_s_1d(ji,1)        = e_s_1d(ji,1)        * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
+                  
+                  ! new concentration
+                  a_i_1d(ji) = a_i_1d(ji) + zda
+                  
+                  ! ensure that ht_i = 0 where a_i = 0
+                  IF( a_i_1d(ji) == 0._wp ) THEN
+                     ht_i_1d(ji) = 0._wp
+                     ht_s_1d(ji) = 0._wp
+                  ENDIF
+               ENDIF
+            END DO
+            
             CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d    (1:nidx), a_i (:,:,jl) )
             CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d   (1:nidx), ht_i(:,:,jl) )
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8369)
@@ -631,5 +631,5 @@
          IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
             sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean
-               &                            - sm_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean 
+               &                            - sm_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean 
          ENDIF
 
@@ -665,6 +665,8 @@
 
       ! --- ensure that a_i = 0 where ht_i = 0 ---
-      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
-      
+      DO ji = 1, nidx
+         IF( ht_i_1d(ji) == 0._wp )   a_i_1d(ji) = 0._wp
+      END DO
+         
       CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i )
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8369)
@@ -735,8 +735,11 @@
       DO ji = 1, nidx
          zfac        = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
-         t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + &
-            &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac
-      END DO
-      WHERE( ( zh_s < 1.0e-3 ) ) ; t_si_1d(:) = t_su_1d(:) ; END WHERE
+         IF( zh_s(ji) >= 1.e-3 ) THEN
+            t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + &
+               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac
+         ELSE
+            t_si_1d(ji) = t_su_1d(ji)
+         ENDIF
+      END DO
       ! END MV SIMIP 2016
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8368)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8369)
@@ -12,5 +12,5 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE ice, ONLY :   nlay_i, nlay_s
+   USE ice, ONLY :   nlay_i, nlay_s, jpl
 
    IMPLICIT NONE
@@ -19,7 +19,7 @@
    PUBLIC thd_ice_alloc ! Routine called by nemogcm.F90
 
-   !!-----------------------------
-   !! * Share 1D Module variables
-   !!-----------------------------
+   !!----------------------
+   !! * 1D Module variables
+   !!----------------------
    !: In ice thermodynamics, to spare memory, the vectors are folded
    !: from 1D to 2D vectors. The following variables, with ending _1d
@@ -90,5 +90,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qevap_ice_1d  !: <==> the 3D  qevap_ice
    !                                                     ! to reintegrate longwave flux inside the ice thermodynamics
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice
@@ -98,6 +97,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_si_1d       !: <==> the 2D  t_si
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_s
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_1d       !: <==> the 2D  ht_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ib_1d       !: <==> the 2D  a_i_b
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_ib_1d      !: <==> the 2D  ht_i_b
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_1d       !: <==> the 2D  ht_s
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su         !: Surface Conduction flux 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i       !: Bottom  Conduction flux 
@@ -109,4 +110,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_1d       !: Ice bulk salinity [ppt]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d        !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d        !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_i_1d        !:
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d   !: corresponding to the 2D var  t_s
@@ -127,4 +131,21 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sss_1d
 
+   ! 
+   !!----------------------
+   !! * 2D Module variables
+   !!----------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_i_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_i_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_s_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   oa_i_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smv_i_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_i_2d
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ib_2d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_ib_2d
+   
    INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point
 
@@ -141,5 +162,5 @@
       !!---------------------------------------------------------------------!
       INTEGER ::   thd_ice_alloc   ! return value
-      INTEGER ::   ierr(6), ii
+      INTEGER ::   ierr(7), ii
       !!---------------------------------------------------------------------!
       ierr(:) = 0
@@ -164,5 +185,5 @@
          &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,              &
          &      wfx_lam_1d(jpij)  , dqns_ice_1d(jpij) , evap_ice_1d (jpij),                     &
-         &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &  
+         &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &  
          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
          &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
@@ -170,9 +191,9 @@
       !
       ii = ii + 1
-      ALLOCATE( t_su_1d   (jpij) , t_si_1d  (jpij)  , a_i_1d   (jpij) , ht_i_1d  (jpij) ,    &
-         &      ht_s_1d   (jpij) , fc_su     (jpij) , fc_bo_i  (jpij) ,                      &    
-         &      dh_s_tot  (jpij) , dh_i_surf (jpij) , dh_i_sub (jpij) ,                      &    
-         &      dh_i_bott (jpij) , dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    &
-         &      STAT=ierr(ii) )
+      ALLOCATE( t_su_1d  (jpij) , t_si_1d   (jpij) , a_i_1d  (jpij) , a_ib_1d(jpij) ,                  &
+         &      ht_i_1d  (jpij) , ht_ib_1d  (jpij) , ht_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &    
+         &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &    
+         &      dh_i_bott(jpij) , dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new(jpij) , &
+         &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , STAT=ierr(ii) )
       !
       ii = ii + 1
@@ -186,4 +207,10 @@
       ii = ii + 1
       ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) )
+      !
+      ii = ii + 1
+      ALLOCATE( a_i_2d(jpij,jpl) , a_ib_2d(jpij,jpl) , ht_i_2d(jpij,jpl) , ht_ib_2d(jpij,jpl) , &
+         &      v_i_2d(jpij,jpl) ,v_s_2d(jpij,jpl) ,oa_i_2d(jpij,jpl) ,smv_i_2d(jpij,jpl) ,  &
+         &      a_ip_2d(jpij,jpl) ,v_ip_2d(jpij,jpl) ,t_su_2d(jpij,jpl) ,  &
+         &      STAT=ierr(ii) )
 
       thd_ice_alloc = MAXVAL( ierr(:) )
