Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 5053)
@@ -51,5 +51,4 @@
    REAL(wp), PUBLIC ::   ahi0        !: sea-ice hor. eddy diffusivity coeff. (m2/s)
    REAL(wp), PUBLIC ::   alphaevp    !: coefficient for the solution of EVP int. stresses
-   REAL(wp), PUBLIC ::   hminrhg = 0.001_wp    !: clem : ice volume (a*h in m) below which ice velocity is set to ocean velocity
 
    REAL(wp), PUBLIC ::   usecc2                !:  = 1.0 / ( ecc * ecc )
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 5053)
@@ -228,5 +228,5 @@
          &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
          &                c_rhg, etamn, creepl, ecc, ahi0,                  &
-         &                nevp, telast, alphaevp, hminrhg
+         &                nevp, telast, alphaevp
       !!-------------------------------------------------------------------
                     
@@ -262,5 +262,4 @@
          WRITE(numout,*) '       timescale for elastic waves telast = ', telast
          WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp
-         WRITE(numout,*) '       min ice thickness for rheology calculations     hminrhg = ', hminrhg
       ENDIF
       !
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5053)
@@ -178,5 +178,4 @@
    REAL(wp), PUBLIC ::   ahi0             !: sea-ice hor. eddy diffusivity coeff. (m2/s)
    REAL(wp), PUBLIC ::   relast           !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
-   REAL(wp), PUBLIC ::   hminrhg          !: ice volume (a*h, in m) below which ice velocity is set to ocean velocity
 
    !                                     !!** ice-salinity namelist (namicesal) **
@@ -300,7 +299,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
 
-   ! lateral melting (mono-category)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_melt !: net melting (m)
-
    !!--------------------------------------------------------------------------
    !! * Ice global state variables
@@ -370,4 +366,7 @@
    !!--------------------------------------------------------------------------
    !                                                  !!: ** Namelist namicerun read in sbc_lim_init **
+!clem   INTEGER               , PUBLIC ::   jpl             !: number of ice  categories 
+!clem   INTEGER               , PUBLIC ::   nlay_i          !: number of ice  layers 
+!clem   INTEGER               , PUBLIC ::   nlay_s          !: number of snow layers 
    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
@@ -439,5 +438,5 @@
          &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,  qlead  (jpi,jpj) ,     &
          &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) , &
-         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) , dh_i_melt(jpi,jpj,jpl) , &
+         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl),    &
          &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,   &
Index: anches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90	(revision 5052)
+++ 	(revision )
@@ -1,188 +1,0 @@
-MODULE limcat_1D
-   !!======================================================================
-   !!                     ***  MODULE  limcat_1D  ***
-   !!  Used for LIM3 to convert cell averages of ice thickness, snow thickness
-   !!  and ice cover into a prescribed distribution over the cell. 
-   !!  (Example of application: BDY forcings when input are cell averaged)  
-   !!======================================================================
-   !! History :   -   ! Original code from M. Vancoppenolle (?)
-   !!                 ! 2011-12 (C. Rousset) rewritten for clarity
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                    LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_cat_1D      :  main subroutine
-   !!----------------------------------------------------------------------
-   !! Modules used
-   USE phycst
-   USE oce             ! dynamics and tracers variables
-   USE dom_oce
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE par_ice         ! ice parameters
-   USE ice             ! ice variables
-   USE eosbn2          ! equation of state
-   USE in_out_manager
-   USE dom_ice
-   USE ice
-   USE lbclnk
-   USE timing          ! Timing
-   USE wrk_nemo
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! Accessibility
-   PUBLIC lim_cat_1D      
-
-CONTAINS
-
-   SUBROUTINE lim_cat_1D( zhti, zhts, zai, zht_i, zht_s, za_i )
-      !! Local variables
-      INTEGER  :: ji, jk, jl             ! dummy loop indices
-      INTEGER  :: ijpij, i_fill, jl0  
-      REAL(wp) :: zarg, zV, zconv, zdh
-      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
-      INTEGER , POINTER, DIMENSION(:)         ::   itest
- 
-      IF( nn_timing == 1 )  CALL timing_start('limcat_1D')
-
-      CALL wrk_alloc( 4, itest )
-      !--------------------------------------------------------------------
-      ! initialisation of variables
-      !--------------------------------------------------------------------
-      ijpij = SIZE(zhti,1)
-      zht_i(1:ijpij,1:jpl) = 0._wp
-      zht_s(1:ijpij,1:jpl) = 0._wp
-      za_i (1:ijpij,1:jpl) = 0._wp
-
-      !------------------------------------------------------------------------------------
-      ! Distribute ice concentration and thickness into the categories
-      !------------------------------------------------------------------------------------
-      ! Method: we first try to fill the jpl ice categories bounded by thicknesses
-      !         hmax(0:jpl) with a gaussian distribution, and check whether the distribution
-      !         fulfills volume and area conservation, positivity and ice categories bounds.
-      !         In other words, if ice input is too thin, the last category (jpl)
-      !         cannot be filled, so we try to fill jpl-1 categories...
-      !         And so forth iteratively until the number of categories filled
-      !         fulfills ice volume concervation between input and output (itest=4) 
-      !--------------------------------------------------------------------------------------
-
-      ! ----------------------------------------
-      ! distribution over the jpl ice categories
-      ! ----------------------------------------
-      DO ji = 1, ijpij
-         
-         IF( zhti(ji) > 0._wp ) THEN
-
-         ! initialisation of tests
-         itest(:)  = 0
-         
-         i_fill = jpl + 1                                             !====================================
-         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories  
-            ! iteration                                               !====================================
-            i_fill = i_fill - 1
-            
-            ! initialisation of ice variables for each try
-            zht_i(ji,1:jpl) = 0._wp
-            za_i (ji,1:jpl) = 0._wp
-            
-            ! *** case very thin ice: fill only category 1
-            IF ( i_fill == 1 ) THEN
-               zht_i(ji,1) = zhti(ji)
-               za_i (ji,1) = zai (ji)
-
-            ! *** case ice is thicker: fill categories >1
-            ELSE
-
-               ! Fill ice thicknesses except the last one (i_fill) by hmean 
-               DO jl = 1, i_fill - 1
-                  zht_i(ji,jl) = hi_mean(jl)
-               END DO
-               
-               ! find which category (jl0) the input ice thickness falls into
-               jl0 = i_fill
-               DO jl = 1, i_fill
-                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
-                     jl0 = jl
-		     CYCLE
-                  ENDIF
-               END DO
-               
-               ! Concentrations in the (i_fill-1) categories 
-               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
-               DO jl = 1, i_fill - 1
-                  IF ( jl == jl0 ) CYCLE
-                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
-                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
-               END DO
-               
-               ! Concentration in the last (i_fill) category
-               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
-               
-               ! Ice thickness in the last (i_fill) category
-               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
-               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
-               
-            ENDIF ! case ice is thick or thin
-            
-            !---------------------
-            ! Compatibility tests
-            !--------------------- 
-            ! Test 1: area conservation
-            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
-            IF ( zconv < epsi06 ) itest(1) = 1
-            
-            ! Test 2: volume conservation
-            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
-            IF ( zconv < epsi06 ) itest(2) = 1
-            
-            ! Test 3: thickness of the last category is in-bounds ?
-            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
-            
-            ! Test 4: positivity of ice concentrations
-            itest(4) = 1
-            DO jl = 1, i_fill
-               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
-            END DO            
-                                                           !============================
-         END DO                                            ! end iteration on categories
-                                                           !============================
-         ENDIF ! if zhti > 0
-      END DO ! i loop
-
-      ! ------------------------------------------------
-      ! Adding Snow in each category where za_i is not 0
-      ! ------------------------------------------------ 
-      DO jl = 1, jpl
-         DO ji = 1, ijpij
-            IF( za_i(ji,jl) > 0._wp ) THEN
-               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
-               ! In case snow load is in excess that would lead to transformation from snow to ice
-               ! Then, transfer the snow excess into the ice (different from limthd_dh)
-               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
-               ! recompute ht_i, ht_s avoiding out of bounds values
-               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
-               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
-            ENDIF
-         ENDDO
-      ENDDO
-
-      CALL wrk_dealloc( 4, itest )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D')
-     
-    END SUBROUTINE lim_cat_1D
-    
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :         Empty module          NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_cat_1D          ! Empty routine
-   END SUBROUTINE lim_cat_1D
-#endif
-
-   !!======================================================================
-END MODULE limcat_1D
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5053)
@@ -31,4 +31,5 @@
    USE timing          ! Timing
    USE limcons        ! conservation tests
+   USE limvar
 
    IMPLICIT NONE
@@ -75,4 +76,6 @@
       CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
       CALL wrk_alloc( jpj, zswitch, zmsk )
+
+      CALL lim_var_agg(1)             ! aggregate ice categories
 
       IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
@@ -241,7 +244,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ cw, pstar,   &
-         &                c_rhg, creepl, ecc, ahi0,     &
-         &                nevp, relast, hminrhg
+      NAMELIST/namicedyn/ cw, pstar, c_rhg, creepl, ecc, ahi0, nevp, relast
       !!-------------------------------------------------------------------
 
@@ -267,5 +268,4 @@
          WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
          WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast
-         WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg
       ENDIF
       !
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5053)
@@ -375,4 +375,5 @@
       CALL lim_var_glo2eqv
       CALL lim_var_zapsmall
+      CALL lim_var_agg( 1 ) 
 
 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5053)
@@ -13,5 +13,4 @@
    !!   'key_lim3' :                                   LIM3 sea-ice model
    !!----------------------------------------------------------------------
-   !!   lim_itd_th       : thermodynamics of ice thickness distribution
    !!   lim_itd_th_rem   :
    !!   lim_itd_th_reb   :
@@ -26,5 +25,4 @@
    USE ice              ! LIM-3 variables
    USE par_ice          ! LIM-3 parameters
-   USE limthd_lac       ! LIM-3 lateral accretion
    USE limvar           ! LIM-3 variables
    USE limcons          ! LIM-3 conservation
@@ -34,11 +32,9 @@
    USE wrk_nemo         ! work arrays
    USE lib_fortran      ! to use key_nosignedzero
-   USE timing          ! Timing
-   USE limcons        ! conservation tests
+   USE limcons          ! conservation tests
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_itd_th         ! called by ice_stp
    PUBLIC   lim_itd_th_rem
    PUBLIC   lim_itd_th_reb
@@ -52,87 +48,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE lim_itd_th( kt )
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_itd_th ***
-      !!
-      !! ** Purpose :   computes the thermodynamics of ice thickness distribution
-      !!
-      !! ** Method  :
-      !!------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! time step index
-      !
-      INTEGER ::   ji, jj, jk, jl   ! dummy loop index         
-      !
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
-      !!------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('limitd_th')
-
-      ! conservation test
-      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_itd_th  : Thermodynamics of the ice thickness distribution'
-         WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-
-      !------------------------------------------------------------------------------|
-      !  1) Transport of ice between thickness categories.                           |
-      !------------------------------------------------------------------------------|
-      ! Given thermodynamic growth rates, transport ice between
-      ! thickness categories.
-      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt )
-      !
-      CALL lim_var_glo2eqv    ! only for info
-      CALL lim_var_agg(1)
-
-      !------------------------------------------------------------------------------|
-      !  3) Add frazil ice growing in leads.
-      !------------------------------------------------------------------------------|
-      CALL lim_thd_lac
-      CALL lim_var_glo2eqv    ! only for info
-
-      ! MV: Could put lateral melting here, would be better I think ???
-
-     
-      IF(ln_ctl) THEN   ! Control print
-         CALL prt_ctl_info(' ')
-         CALL prt_ctl_info(' - Cell values : ')
-         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :')
-         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :')
-         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :')
-         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :')
-         DO jl = 1, jpl
-            CALL prt_ctl_info(' ')
-            CALL prt_ctl_info(' - Category : ', ivar1=jl)
-            CALL prt_ctl_info('   ~~~~~~~~~~')
-            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ')
-            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ')
-            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ')
-            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ')
-            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ')
-            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ')
-            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ')
-            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ')
-            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ')
-            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ')
-            DO jk = 1, nlay_i
-               CALL prt_ctl_info(' ')
-               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
-               CALL prt_ctl_info('   ~~~~~~~')
-               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ')
-               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ')
-            END DO
-         END DO
-      ENDIF
-      !
-      ! conservation test
-      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      !
-     IF( nn_timing == 1 )  CALL timing_stop('limitd_th')
-   END SUBROUTINE lim_itd_th
-   !
 
    SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
@@ -156,5 +69,5 @@
       REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
       REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
-      REAL(wp) ::   zx3,             zareamin          !   -      -
+      REAL(wp) ::   zx3        
       CHARACTER (len = 15) :: fieldid
 
@@ -191,6 +104,4 @@
       CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
 
-      zareamin = epsi10   !minimum area in thickness categories tolerated by the conceptors of the model
-
       !!----------------------------------------------------------------------------------------------
       !! 0) Conservation checkand changes in each ice category
@@ -242,5 +153,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF ( at_i(ji,jj) .gt. zareamin ) THEN
+            IF ( at_i(ji,jj) > epsi10 ) THEN
                nbrem         = nbrem + 1
                nind_i(nbrem) = ji
@@ -1014,8 +925,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE lim_itd_th           ! Empty routines
-   END SUBROUTINE lim_itd_th
-   SUBROUTINE lim_itd_th_ini
-   END SUBROUTINE lim_itd_th_ini
    SUBROUTINE lim_itd_th_rem
    END SUBROUTINE lim_itd_th_rem
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5053)
@@ -152,4 +152,5 @@
 
       REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter
+      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
       !!-------------------------------------------------------------------
 
@@ -161,14 +162,14 @@
 #if  defined key_lim2 && ! defined key_lim2_vp
 # if defined key_agrif
-     USE ice_2, vt_s => hsnm
-     USE ice_2, vt_i => hicm
+      USE ice_2, vt_s => hsnm
+      USE ice_2, vt_i => hicm
 # else
-     vt_s => hsnm
-     vt_i => hicm
+      vt_s => hsnm
+      vt_i => hicm
 # endif
-     at_i(:,:) = 1. - frld(:,:)
+      at_i(:,:) = 1. - frld(:,:)
 #endif
 #if defined key_agrif && defined key_lim2 
-    CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
+      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
 #endif
       !
@@ -189,7 +190,5 @@
 #endif
 
-!CDIR NOVERRCHK
       DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
-!CDIR NOVERRCHK
          DO ji = 1 , jpi
 #if defined key_lim3
@@ -206,20 +205,13 @@
       ! Ice strength on grid cell corners (zpreshc)
       ! needed for calculation of shear stress 
-!CDIR NOVERRCHK
       DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
          DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
-            zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
-               &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
-               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &
-               &              tms(ji,jj)     * wght(ji+1,jj+1,1,1)
-            zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
-               &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
-               &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
-               &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   &
+            zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
+               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + tms(ji,jj)   * wght(ji+1,jj+1,1,1)
+            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
+               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     &
                &             ) / MAX( zstms, zepsi )
          END DO
       END DO
-
       CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
       !
@@ -242,15 +234,15 @@
 
       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
-          !                                            
-          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
-          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
+         !                                            
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
          zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
-          !
-          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
-          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         !
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
-          !
+         !
          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
-          !
+         !
       ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
          zpice(:,:) = ssh_m(:,:)
@@ -363,40 +355,26 @@
                !
                !
-               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
-                  &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
-                  &             +e1v(ji,jj)*v_ice(ji,jj)                      &
-                  &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
-                  &             )                                             &
-                  &            / area(ji,jj)
-
-               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
-                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
-                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
-                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
-                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
-                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
-                  &         )                                              &
-                  &        / area(ji,jj)
+               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &            ) / area(ji,jj)
+
+               zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &         - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &         ) / area(ji,jj)
 
                !
-               zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                &
-                  &            -u_ice(ji,jj)/e1u(ji,jj)                    &
-                  &           )*e1f(ji,jj)*e1f(ji,jj)                      &
-                  &          +( v_ice(ji+1,jj)/e2v(ji+1,jj)                &
-                  &            -v_ice(ji,jj)/e2v(ji,jj)                    &
-                  &           )*e2f(ji,jj)*e2f(ji,jj)                      &
-                  &         )                                              &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
-                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
-                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
-
-
-               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
-                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
-                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
-
-               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
-                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
-                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+                  &         + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+                  &         ) / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2._wp - tmf(ji,jj) )   &
+                  &         * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1)
+
+
+               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji,jj)   + v_ice(ji,jj-1)   ) * e1t(ji+1,jj)   &
+                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji,jj) )   &
+                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 
+
+               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj)   + u_ice(ji-1,jj)   ) * e2t(ji,jj+1)   &
+                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj) )   &
+                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)
 
             END DO
@@ -404,18 +382,13 @@
          CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond.
 
-!CDIR NOVERRCHK
          DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1
 
                !- Calculate Delta at centre of grid cells
-               zdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
-                  &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          &
-                  &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          &
-                  &          - e1v(ji, jj-1) * u_ice2(ji,jj-1)          &
-                  &          )                                          &
-                  &         / area(ji,jj)
-
-               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
+               zdst      = ( e2u(ji, jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
+                  &        + e1v(ji, jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
+                  &        ) / area(ji,jj)
+
+               delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) * usecc2 )  
                delta_i(ji,jj) = delta + creepl
                !-Calculate stress tensor components zs1 and zs2 
@@ -432,7 +405,5 @@
          CALL lbc_lnk( zs2(:,:), 'T', 1. )
 
-!CDIR NOVERRCHK
          DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1
                !- Calculate Delta on corners
@@ -490,7 +461,5 @@
          IF (MOD(jter,2).eq.0) THEN 
 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
@@ -520,7 +489,5 @@
 #endif         
 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
 
@@ -551,7 +518,5 @@
 
          ELSE 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
@@ -581,7 +546,5 @@
 #endif         
 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
@@ -628,12 +591,9 @@
       ! 4) Prevent ice velocities when the ice is thin
       !------------------------------------------------------------------------------!
-      ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the
-      ! ocean velocity, 
-      ! This prevents high velocity when ice is thin
-!CDIR NOVERRCHK
+      ! If the ice volume is below zvmin then ice velocity should equal the
+      ! ocean velocity. This prevents high velocity when ice is thin
       DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1
-            IF ( vt_i(ji,jj) .LE. hminrhg ) THEN
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
                u_ice(ji,jj) = u_oce(ji,jj)
                v_ice(ji,jj) = v_oce(ji,jj)
@@ -655,12 +615,12 @@
       DO jj = k_j1+1, k_jpj-1 
          DO ji = fs_2, fs_jpim1
-            IF ( vt_i(ji,jj) .LE. hminrhg ) THEN
-               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
-                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
-                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
-
-               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
-                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
-                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
+               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     &
+                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
+                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)
+
+               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
+                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)
             ENDIF 
          END DO
@@ -671,11 +631,9 @@
 
       ! Recompute delta, shear and div, inputs for mechanical redistribution 
-!CDIR NOVERRCHK
       DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
          DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
             !- divu_i(:,:), zdt(:,:): divergence and tension at centre 
             !- zds(:,:): shear on northeast corner of grid cells
-            IF ( vt_i(ji,jj) .LE. hminrhg ) THEN
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
 
                divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
@@ -722,10 +680,8 @@
       ! 5) Store stress tensor and its invariants
       !------------------------------------------------------------------------------!
-      !
       ! * Invariants of the stress tensor are required for limitd_me
       !   (accelerates convergence and improves stability)
       DO jj = k_j1+1, k_jpj-1
          DO ji = fs_2, fs_jpim1
-            ! begin TECLIM change 
             zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
                &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &   
@@ -733,5 +689,4 @@
                &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 
             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
-            ! end TECLIM change
          END DO
       END DO
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 5053)
@@ -27,4 +27,5 @@
    USE wrk_nemo       ! work arrays
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE limctl
 
    IMPLICIT NONE
@@ -86,4 +87,5 @@
       ENDIF
       !
+      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
    END SUBROUTINE lim_rst_opn
 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5053)
@@ -43,4 +43,5 @@
    USE iom
    USE domvvl           ! Variable volume
+   USE limctl
 
    IMPLICIT NONE
@@ -225,4 +226,5 @@
       ENDIF
 
+      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
 
       IF(ln_ctl) THEN
@@ -270,7 +272,5 @@
       !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
-!CDIR NOVERRCHK
          DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1
                !                                               ! 2*(U_ice-U_oce) at T-point
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5053)
@@ -34,4 +34,6 @@
    USE limthd_sal     ! LIM: thermodynamics, ice salinity
    USE limthd_ent     ! LIM: thermodynamics, ice enthalpy redistribution
+   USE limthd_lac     ! LIM-3 lateral accretion
+   USE limitd_th      ! remapping thickness distribution
    USE limtab         ! LIM: 1D <==> 2D transformation
    USE limvar         ! LIM: sea-ice variables
@@ -44,4 +46,5 @@
    USE timing         ! Timing
    USE limcons        ! conservation tests
+   USE limctl
 
    IMPLICIT NONE
@@ -80,16 +83,13 @@
       !! ** References : 
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! number of iteration
+      INTEGER, INTENT(in) :: kt    ! number of iteration
       !!
       INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
       INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations
-      INTEGER  :: nbplm            ! nb of icy pts for lateral melting calculations (mono-cat)
       INTEGER  :: ii, ij           ! temporary dummy loop index
-      REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04)
-      REAL(wp) :: zch        = 0.0057_wp    ! heat transfer coefficient
-      REAL(wp) :: zareamin 
       REAL(wp) :: zfric_u, zqld, zqfr
-      !
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      REAL(wp), PARAMETER :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04)
+      REAL(wp), PARAMETER :: zch        = 0.0057_wp    ! heat transfer coefficient
       !
       REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns
@@ -106,5 +106,4 @@
       !------------------------------------------------------------------------!
       ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice
-
 
       !--------------------
@@ -162,9 +161,7 @@
       ENDIF
 
-!CDIR NOVERRCHK
       DO jj = 1, jpj
-!CDIR NOVERRCHK
          DO ji = 1, jpi
-            rswitch          = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice
+            rswitch  = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice
             !
             !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
@@ -260,9 +257,8 @@
          ENDIF
 
-         zareamin = epsi10
          nbpb = 0
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN     
+               IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
                   nbpb      = nbpb  + 1
                   npb(nbpb) = (jj - 1) * jpi + ji
@@ -442,8 +438,4 @@
               CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj )
               CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
-
-!clem            IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN
-!clem              CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj )
-!clem            ENDIF
             !
               CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
@@ -460,5 +452,5 @@
 
       !------------------------
-      ! 5.1) Ice heat content              
+      ! Ice heat content              
       !------------------------
       ! Enthalpies are global variables we have to readjust the units (heat content in Joules)
@@ -470,5 +462,5 @@
 
       !------------------------
-      ! 5.2) Snow heat content              
+      ! Snow heat content              
       !------------------------
       ! Enthalpies are global variables we have to readjust the units (heat content in Joules)
@@ -478,7 +470,12 @@
          END DO
       END DO
+ 
+      !------------------------
+      ! Ice natural aging              
+      !------------------------
+      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday
 
       !----------------------------------
-      ! 5.3) Change thickness to volume
+      ! Change thickness to volume
       !----------------------------------
       CALL lim_var_eqv2glo
@@ -487,4 +484,6 @@
       ! 5.4) Diagnostic thermodynamic growth rates
       !--------------------------------------------
+      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
+
       IF(ln_ctl) THEN            ! Control print
          CALL prt_ctl_info(' ')
@@ -522,8 +521,55 @@
       CALL wrk_dealloc( jpi, jpj, zqsr, zqns )
 
+      !------------------------------------------------------------------------------|
+      !  1) Transport of ice between thickness categories.                           |
+      !------------------------------------------------------------------------------|
+      ! Given thermodynamic growth rates, transport ice between
+      ! thickness categories.
+      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt )
+      !
+      CALL lim_var_glo2eqv    ! only for info
+      CALL lim_var_agg(1)
+
+      !------------------------------------------------------------------------------|
+      !  3) Add frazil ice growing in leads.
+      !------------------------------------------------------------------------------|
+      CALL lim_thd_lac
+      CALL lim_var_glo2eqv    ! only for info
+      
+      IF(ln_ctl) THEN   ! Control print
+         CALL prt_ctl_info(' ')
+         CALL prt_ctl_info(' - Cell values : ')
+         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
+         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :')
+         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :')
+         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :')
+         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :')
+         DO jl = 1, jpl
+            CALL prt_ctl_info(' ')
+            CALL prt_ctl_info(' - Category : ', ivar1=jl)
+            CALL prt_ctl_info('   ~~~~~~~~~~')
+            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ')
+            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ')
+            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ')
+            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ')
+            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ')
+            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ')
+            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ')
+            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ')
+            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ')
+            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ')
+            DO jk = 1, nlay_i
+               CALL prt_ctl_info(' ')
+               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
+               CALL prt_ctl_info('   ~~~~~~~')
+               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ')
+               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ')
+            END DO
+         END DO
+      ENDIF
       !
       ! conservation test
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      !
+
       IF( nn_timing == 1 )  CALL timing_stop('limthd')
 
@@ -571,10 +617,12 @@
       INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
       INTEGER             ::   ji            ! dummy loop indices
-
-      WRITE(numout,*) ' Lateral melting ON '
+      REAL(wp)            ::   zhi_bef       ! ice thickness before thermo
+      REAL(wp)            ::   zdh_mel       ! net melting
+
       DO ji = kideb, kiut
-         IF( ht_i_1d(ji) > epsi10 .AND. dh_i_melt_1d(ji) < 0._wp ) THEN     
-            a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) ) ) 
-         END IF
+         zhi_bef = ht_i_1d(ji) - ( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
+         zdh_mel = MIN( 0._wp,     dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
+         IF( zdh_mel < 0._wp )   &
+            &   a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) ) 
       END DO
       at_i_1d(:) = a_i_1d(:)
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5053)
@@ -689,15 +689,4 @@
       END DO
       
-      !
-      !-------------------------------------------------------------
-      ! Net melting (input for thin ice melting, mono-category case)
-      !-------------------------------------------------------------
-      !
-      IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN
-         DO ji = kideb, kiut
-            dh_i_melt_1d(ji) = MIN( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji), 0._wp )
-         END DO
-      ENDIF
-
       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema )
       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5053)
@@ -22,5 +22,5 @@
    USE limadv         ! ice advection
    USE limhdf         ! ice horizontal diffusion
-   USE limvar         ! clem for ice thickness correction
+   USE limvar         ! 
    !
    USE in_out_manager ! I/O manager
@@ -32,4 +32,5 @@
    USE timing         ! Timing
    USE limcons        ! conservation tests
+   USE limctl         ! control prints
 
    IMPLICIT NONE
@@ -455,4 +456,9 @@
       ENDIF
 
+      CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
+
+      ! -------------------------------------------------
+      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
+      ! -------------------------------------------------
       IF(ln_ctl) THEN   ! Control print
          CALL prt_ctl_info(' ')
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5053)
@@ -59,5 +59,5 @@
 CONTAINS
 
-   SUBROUTINE lim_update1
+   SUBROUTINE lim_update1( kt )
       !!-------------------------------------------------------------------
       !!               ***  ROUTINE lim_update1  ***
@@ -67,4 +67,5 @@
       !!                
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt    ! number of iteration
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       INTEGER  ::   i_ice_switch
@@ -76,4 +77,9 @@
 
       IF( ln_limdyn ) THEN 
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*) ' lim_update1 ' 
+         WRITE(numout,*) ' ~~~~~~~~~~~ '
+      ENDIF
 
       ! conservation test
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5053)
@@ -41,4 +41,5 @@
    USE timing          ! Timing
    USE limcons        ! conservation tests
+   USE limctl
 
    IMPLICIT NONE
@@ -56,5 +57,5 @@
 CONTAINS
 
-   SUBROUTINE lim_update2
+   SUBROUTINE lim_update2( kt )
       !!-------------------------------------------------------------------
       !!               ***  ROUTINE lim_update2  ***
@@ -64,4 +65,5 @@
       !!
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt    ! number of iteration
       INTEGER  ::   ji, jj, jk, jl    ! dummy loop indices
       INTEGER  ::   i_ice_switch
@@ -71,4 +73,9 @@
       !!-------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limupdate2')
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*) ' lim_update2 '
+         WRITE(numout,*) ' ~~~~~~~~~~~ '
+      ENDIF
 
       ! conservation test
@@ -208,4 +215,6 @@
       ! conservation test
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print
 
       IF(ln_ctl) THEN   ! Control print
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5053)
@@ -30,16 +30,9 @@
    !!======================================================================
    !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
    !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
    !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_var_agg       : 
-   !!   lim_var_glo2eqv   :
-   !!   lim_var_eqv2glo   :
-   !!   lim_var_salprof   : 
-   !!   lim_var_salprof1d :
-   !!   lim_var_bv        :
    !!----------------------------------------------------------------------
    USE par_oce        ! ocean parameters
@@ -58,15 +51,16 @@
    PRIVATE
 
-   PUBLIC   lim_var_agg          !
-   PUBLIC   lim_var_glo2eqv      !
-   PUBLIC   lim_var_eqv2glo      !
-   PUBLIC   lim_var_salprof      !
-   PUBLIC   lim_var_icetm        !
-   PUBLIC   lim_var_bv           !
-   PUBLIC   lim_var_salprof1d    !
+   PUBLIC   lim_var_agg          
+   PUBLIC   lim_var_glo2eqv      
+   PUBLIC   lim_var_eqv2glo      
+   PUBLIC   lim_var_salprof      
+   PUBLIC   lim_var_icetm        
+   PUBLIC   lim_var_bv           
+   PUBLIC   lim_var_salprof1d    
    PUBLIC   lim_var_zapsmall
+   PUBLIC   lim_var_itd
 
    !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -130,5 +124,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content
+                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content
                   rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity
@@ -192,17 +186,13 @@
       ! Ice temperatures
       !-------------------
-!CDIR NOVERRCHK
-      DO jl = 1, jpl
-!CDIR NOVERRCHK
+      DO jl = 1, jpl
          DO jk = 1, nlay_i
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
+            DO jj = 1, jpj
                DO ji = 1, jpi
                   !                                                              ! Energy of melting q(S,T) [J.m-3]
                   rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes
                   zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
-                  zq_i    = zq_i * unit_fac                             !convert units
-                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature
+                  zq_i    = zq_i * unit_fac                             ! convert units
+                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt              ! Ice layer melt temperature
                   !
                   zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
@@ -289,5 +279,5 @@
       !!              alpha-weighted linear combination of s_inf and s_zero
       !!
-      !! ** References : Vancoppenolle et al., 2007 (in preparation)
+      !! ** References : Vancoppenolle et al., 2007
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
@@ -368,5 +358,4 @@
          !
          DO jl = 1, jpl
-!CDIR NOVERRCHK
             DO jk = 1, nlay_i
                zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
@@ -451,5 +440,5 @@
       !
       INTEGER  ::   ji, jk    ! dummy loop indices
-      INTEGER  ::   ii, ij  ! local integers
+      INTEGER  ::   ii, ij    ! local integers
       REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars
       REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      -
@@ -483,7 +472,5 @@
          dummy_fac2 = 1._wp / REAL(nlay_i,wp)
 
-!CDIR NOVERRCHK
          DO jk = 1, nlay_i
-!CDIR NOVERRCHK
             DO ji = kideb, kiut
                ii =  MOD( npb(ji) - 1 , jpi ) + 1
@@ -502,8 +489,8 @@
                ! weighting the profile
                s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
-            END DO ! ji
-         END DO ! jk
-
-      ENDIF ! num_sal
+            END DO 
+         END DO 
+
+      ENDIF 
 
       !-------------------------------------------------------
@@ -515,5 +502,4 @@
          sm_i_1d(:) = 2.30_wp
          !
-!CDIR NOVERRCHK
          DO jk = 1, nlay_i
             zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
@@ -603,4 +589,160 @@
       !
    END SUBROUTINE lim_var_zapsmall
+
+   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE lim_var_itd   ***
+      !!
+      !! ** Purpose :  converting 1-cat ice to multiple ice categories
+      !!
+      !!                  ice thickness distribution follows a gaussian law
+      !!               around the concentration of the most likely ice thickness
+      !!                           (similar as limistate.F90)
+      !!
+      !! ** Method:   Iterative procedure
+      !!                
+      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
+      !!
+      !!               2) Check whether the distribution conserves area and volume, positivity and
+      !!                  category boundaries
+      !!              
+      !!               3) If not (input ice is too thin), the last category is empty and
+      !!                  the number of categories is reduced (jpl-1)
+      !!
+      !!               4) Iterate until ok (SUM(itest(:) = 4)
+      !!
+      !! ** Arguments : zhti: 1-cat ice thickness
+      !!                zhts: 1-cat snow depth
+      !!                zai : 1-cat ice concentration
+      !!
+      !! ** Output    : jpl-cat 
+      !!
+      !!  (Example of application: BDY forcings when input are cell averaged)  
+      !!
+      !!-------------------------------------------------------------------
+      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
+      !!                    2014    (C. Rousset)        Rewriting
+      !!-------------------------------------------------------------------
+      !! Local variables
+      INTEGER  :: ji, jk, jl             ! dummy loop indices
+      INTEGER  :: ijpij, i_fill, jl0  
+      REAL(wp) :: zarg, zV, zconv, zdh
+      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
+      INTEGER , POINTER, DIMENSION(:)         ::   itest
+ 
+      CALL wrk_alloc( 4, itest )
+      !--------------------------------------------------------------------
+      ! initialisation of variables
+      !--------------------------------------------------------------------
+      ijpij = SIZE(zhti,1)
+      zht_i(1:ijpij,1:jpl) = 0._wp
+      zht_s(1:ijpij,1:jpl) = 0._wp
+      za_i (1:ijpij,1:jpl) = 0._wp
+
+      ! ----------------------------------------
+      ! distribution over the jpl ice categories
+      ! ----------------------------------------
+      DO ji = 1, ijpij
+         
+         IF( zhti(ji) > 0._wp ) THEN
+
+         ! initialisation of tests
+         itest(:)  = 0
+         
+         i_fill = jpl + 1                                             !====================================
+         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories  
+            ! iteration                                               !====================================
+            i_fill = i_fill - 1
+            
+            ! initialisation of ice variables for each try
+            zht_i(ji,1:jpl) = 0._wp
+            za_i (ji,1:jpl) = 0._wp
+            
+            ! *** case very thin ice: fill only category 1
+            IF ( i_fill == 1 ) THEN
+               zht_i(ji,1) = zhti(ji)
+               za_i (ji,1) = zai (ji)
+
+            ! *** case ice is thicker: fill categories >1
+            ELSE
+
+               ! Fill ice thicknesses except the last one (i_fill) by hmean 
+               DO jl = 1, i_fill - 1
+                  zht_i(ji,jl) = hi_mean(jl)
+               END DO
+               
+               ! find which category (jl0) the input ice thickness falls into
+               jl0 = i_fill
+               DO jl = 1, i_fill
+                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
+                     jl0 = jl
+		     CYCLE
+                  ENDIF
+               END DO
+               
+               ! Concentrations in the (i_fill-1) categories 
+               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
+               DO jl = 1, i_fill - 1
+                  IF ( jl == jl0 ) CYCLE
+                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
+                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
+               END DO
+               
+               ! Concentration in the last (i_fill) category
+               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
+               
+               ! Ice thickness in the last (i_fill) category
+               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
+               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
+               
+            ENDIF ! case ice is thick or thin
+            
+            !---------------------
+            ! Compatibility tests
+            !--------------------- 
+            ! Test 1: area conservation
+            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
+            IF ( zconv < epsi06 ) itest(1) = 1
+            
+            ! Test 2: volume conservation
+            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
+            IF ( zconv < epsi06 ) itest(2) = 1
+            
+            ! Test 3: thickness of the last category is in-bounds ?
+            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
+            
+            ! Test 4: positivity of ice concentrations
+            itest(4) = 1
+            DO jl = 1, i_fill
+               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
+            END DO            
+                                                           !============================
+         END DO                                            ! end iteration on categories
+                                                           !============================
+         ENDIF ! if zhti > 0
+      END DO ! i loop
+
+      ! ------------------------------------------------
+      ! Adding Snow in each category where za_i is not 0
+      ! ------------------------------------------------ 
+      DO jl = 1, jpl
+         DO ji = 1, ijpij
+            IF( za_i(ji,jl) > 0._wp ) THEN
+               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
+               ! In case snow load is in excess that would lead to transformation from snow to ice
+               ! Then, transfer the snow excess into the ice (different from limthd_dh)
+               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
+               ! recompute ht_i, ht_s avoiding out of bounds values
+               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
+               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
+            ENDIF
+         ENDDO
+      ENDDO
+
+      CALL wrk_dealloc( 4, itest )
+      !
+    END SUBROUTINE lim_var_itd
+
 
 #else
@@ -623,4 +765,6 @@
    SUBROUTINE lim_var_zapsmall
    END SUBROUTINE lim_var_zapsmall
+   SUBROUTINE lim_var_itd
+   END SUBROUTINE lim_var_itd
 #endif
 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5053)
@@ -110,5 +110,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_melt_1d  !: Net melting [m], for mono-cat lateral melting
    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
@@ -167,5 +166,5 @@
          &      ht_s_1d    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &    
          &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &    
-         &      dh_snowice(jpij) , dh_i_melt_1d(jpij) , &
+         &      dh_snowice(jpij) ,  &
          &      sm_i_1d   (jpij) , s_i_new  (jpij) ,    &
          &      t_s_1d(jpij,nlay_s),                                       &
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 5053)
@@ -35,5 +35,5 @@
    USE par_ice
    USE ice
-   USE limcat_1D          ! redistribute ice input into categories
+   USE limvar          ! redistribute ice input into categories
 #endif
    USE sbcapr
@@ -380,5 +380,5 @@
 #if defined key_lim3
                IF( .NOT. ll_bdylim3 .AND. cn_ice_lim(ib_bdy) /= 'none' .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type)
-     	          CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
+     	          CALL lim_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
                                   & dta_bdy(ib_bdy)%ht_i,     dta_bdy(ib_bdy)%ht_s,     dta_bdy(ib_bdy)%a_i     )
                ENDIF
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 5053)
@@ -60,6 +60,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) :: kt     ! Main time step counter
-      !!
       INTEGER               :: ib_bdy ! Loop index
+
       DO ib_bdy=1, nb_bdy
 
@@ -72,5 +72,6 @@
             CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' )
          END SELECT
-      ENDDO
+
+      END DO
 
    END SUBROUTINE bdy_ice_lim
@@ -247,8 +248,7 @@
             END DO
 
-
-         END DO !jb
+         END DO
  
-         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions
+         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
          CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
          CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5052)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5053)
@@ -39,5 +39,4 @@
    USE limtrp          ! Ice transport
    USE limthd          ! Ice thermodynamics
-   USE limitd_th       ! Thermodynamics on ice thickness distribution 
    USE limitd_me       ! Mechanics on ice thickness distribution
    USE limsbc          ! sea surface boundary condition
@@ -109,6 +108,5 @@
       INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED)
       !!
-      INTEGER  ::   jl      ! dummy loop index
-      REAL(wp) ::   zcoef   ! local scalar
+      INTEGER  ::   jl                 ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled)
@@ -117,10 +115,8 @@
       IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
 
-      !                                        !----------------------!
-      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
-         !                                     !----------------------!
-         !                                           !  Bulk Formulae !
-         !                                           !----------------!
-         !
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only
+         !-----------------------!                                           
+         ! --- Bulk Formulae --- !                                           
+         !-----------------------!
          u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)      ! mean surface ocean current at ice velocity point
          v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)      ! (C-grid dynamics :  U- & V-points as the ocean)
@@ -133,6 +129,7 @@
          CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
 
+         ! CORE and COUPLED bulk formulations
          SELECT CASE( kblk )
-         CASE( jp_core , jp_cpl )   ! CORE and COUPLED bulk formulations
+         CASE( jp_core , jp_cpl )
 
             ! albedo depends on cloud fraction because of non-linear spectral effects
@@ -183,8 +180,7 @@
          END SELECT
          
-         !                                           !----------------------!
-         !                                           ! LIM-3  time-stepping !
-         !                                           !----------------------!
-         ! 
+         !------------------------------!
+         ! --- LIM-3 main time-step --- !
+         !------------------------------!
          numit = numit + nn_fsbc                     ! Ice model time step
          !
@@ -200,9 +196,8 @@
          v_ice_b(:,:)     = v_ice(:,:)
 
-                          CALL sbc_lim_flx0          ! set diag of mass, heat and salt fluxes to 0
+                          CALL sbc_lim_diag0         ! set diag of mass, heat and salt fluxes to 0
 
                           CALL lim_rst_opn( kt )     ! Open Ice restart file
          !
-         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
          ! ----------------------------------------------
          ! ice dynamics and transport (except in 1D case)
@@ -211,10 +206,10 @@
 
                           CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
+
                           CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
-                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
-         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
+
          IF( nn_monocat /= 2 )   &
             &             CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
-                          CALL lim_var_agg( 1 ) 
+
 #if defined key_bdy
                           ! bdy ice thermo 
@@ -225,16 +220,19 @@
          IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print
 #endif
-                          CALL lim_update1
+
+                          CALL lim_update1( kt )
+
          ENDIF
-!                         !- Change old values for new values
-                          u_ice_b(:,:)     = u_ice(:,:)
-                          v_ice_b(:,:)     = v_ice(:,:)
-                          a_i_b  (:,:,:)   = a_i  (:,:,:)
-                          v_s_b  (:,:,:)   = v_s  (:,:,:)
-                          v_i_b  (:,:,:)   = v_i  (:,:,:)
-                          e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
-                          e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
-                          oa_i_b (:,:,:)   = oa_i (:,:,:)
-                          smv_i_b(:,:,:)   = smv_i(:,:,:)
+
+         !- Change old values for new values
+         u_ice_b(:,:)     = u_ice(:,:)
+         v_ice_b(:,:)     = v_ice(:,:)
+         a_i_b  (:,:,:)   = a_i  (:,:,:)
+         v_s_b  (:,:,:)   = v_s  (:,:,:)
+         v_i_b  (:,:,:)   = v_i  (:,:,:)
+         e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
+         e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
+         oa_i_b (:,:,:)   = oa_i (:,:,:)
+         smv_i_b(:,:,:)   = smv_i(:,:,:)
  
          ! ----------------------------------------------
@@ -253,30 +251,21 @@
                              IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
                           &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
-                           ! Latent heat flux is forced to 0 in coupled :
-                           !  it is included in qns (non-solar heat flux)
+                             ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux)
                              qla_ice  (:,:,:) = 0._wp
                              dqla_ice (:,:,:) = 0._wp
                           END SELECT
                           !
-                          CALL lim_var_bv                 ! bulk brine volume (diag)
                           CALL lim_thd( kt )              ! Ice thermodynamics 
-                          zcoef = rdt_ice /rday           !  Ice natural aging
-                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
-         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
-                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion
-                          CALL lim_update2                ! Global variables update
-
-         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print
-         !
-                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes
-         !
-         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
+
+                          CALL lim_update2( kt )          ! Global variables update
+         !
+                          CALL lim_sbc_flx( kt )          ! Update surface ocean mass, heat and salt fluxes
          !
          IF(ln_limdiaout) CALL lim_diahsb                 ! Diagnostics and outputs 
 
-                          CALL lim_wri( 1  )              ! Ice outputs 
+                          CALL lim_wri( 1 )               ! Ice outputs 
 
          IF( kt == nit000 .AND. ln_rstart )   &
-            &             CALL iom_close( numrir )        ! clem: close input ice restart file
+            &             CALL iom_close( numrir )        ! close input ice restart file
          !
          IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file 
@@ -315,5 +304,4 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
       !
-      !
       !                                ! Allocate the ice arrays
       ierr =        ice_alloc        ()      ! ice variables
@@ -343,10 +331,7 @@
       CALL lim_thd_sal_init            ! set ice salinity parameters
       !
-      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep
-      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse
-      !
       CALL lim_msh                     ! ice mesh initialization
       !
-      CALL lim_itd_init                 ! ice thickness distribution initialization
+      CALL lim_itd_init                ! ice thickness distribution initialization
       !
       CALL lim_itd_me_init             ! ice thickness distribution initialization
@@ -356,12 +341,10 @@
          numit = nit000 - 1
          CALL lim_istate
-         CALL lim_var_agg(1)
-         CALL lim_var_glo2eqv
       ELSE                                    ! start from a restart file
          CALL lim_rst_read
          numit = nit000 - 1
-         CALL lim_var_agg(1)
-         CALL lim_var_glo2eqv
       ENDIF
+      CALL lim_var_agg(1)
+      CALL lim_var_glo2eqv
       !
       CALL lim_sbc_init                 ! ice surface boundary condition   
@@ -390,5 +373,8 @@
       !! ** input   :   Namelist namicerun
       !!-------------------------------------------------------------------
-      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout
+!clem      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_out,   &
+!clem         &                ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout
+      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out,   &
+         &                ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!-------------------------------------------------------------------
@@ -408,4 +394,7 @@
          WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
          WRITE(numout,*) ' ~~~~~~'
+!clem         WRITE(numout,*) '   number of ice  categories                               = ', jpl
+!clem         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i
+!clem         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
          WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
          WRITE(numout,*) '   maximum ice concentration                               = ', amax 
@@ -423,4 +412,8 @@
          IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
       ENDIF
+      !
+      ! sea-ice timestep and inverse
+      rdt_ice   = nn_fsbc * rdttra(1)  
+      r1_rdtice = 1._wp / rdt_ice 
       !
    END SUBROUTINE ice_run
@@ -577,5 +570,5 @@
    END SUBROUTINE ice_lim_flx
 
-   SUBROUTINE sbc_lim_flx0
+   SUBROUTINE sbc_lim_diag0
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_lim_flx0  ***
@@ -610,5 +603,5 @@
       afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp
       
-   END SUBROUTINE sbc_lim_flx0
+   END SUBROUTINE sbc_lim_diag0
       
    FUNCTION fice_cell_ave ( ptab )
@@ -623,6 +616,5 @@
       
       DO jl = 1, jpl
-         fice_cell_ave (:,:) = fice_cell_ave (:,:) &
-            &                  + a_i (:,:,jl) * ptab (:,:,jl)
+         fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
       END DO
       
