Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 6801)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 6853)
@@ -31,4 +31,8 @@
    USE iom
 
+!!!clem
+!!   USE diawri
+!!!
+   
    IMPLICIT NONE
    PRIVATE
@@ -84,14 +88,16 @@
       REAL(wp)   :: ztmelts, zdh
       INTEGER    :: i_hemis, i_fill, jl0  
-      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
+      REAL(wp)   :: zarg, zV, zconv, zdv 
       REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator
       REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
       REAL(wp), POINTER, DIMENSION(:,:)   :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini               !data by cattegories to fill
-      !--------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini )
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini                         !data by cattegories to fill
+      INTEGER , POINTER, DIMENSION(:)     :: itest
+      !--------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini )
       CALL wrk_alloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
       CALL wrk_alloc( jpi, jpj,      zswitch )
+      Call wrk_alloc( 4,             itest )
 
       IF(lwp) WRITE(numout,*)
@@ -124,6 +130,6 @@
          DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
             DO ji = 1, jpi
-               IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 
-                  zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
+               IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst * tmask(ji,jj,1) ) THEN
+                  zswitch(ji,jj) = 0._wp                     ! no ice
                ELSE                                                                                   
                   zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
@@ -153,17 +159,17 @@
                DO ji = 1, jpi
                   IF( ff(ji,jj) >= 0._wp ) THEN
-                     zht_i_ini(ji,jj) = rn_hti_ini_n
-                     zht_s_ini(ji,jj) = rn_hts_ini_n
-                     zat_i_ini(ji,jj) = rn_ati_ini_n
-                     zts_u_ini(ji,jj) = rn_tmi_ini_n
-                     zsm_i_ini(ji,jj) = rn_smi_ini_n
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_n
+                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
+                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
+                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
+                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
+                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
                   ELSE
-                     zht_i_ini(ji,jj) = rn_hti_ini_s
-                     zht_s_ini(ji,jj) = rn_hts_ini_s
-                     zat_i_ini(ji,jj) = rn_ati_ini_s
-                     zts_u_ini(ji,jj) = rn_tmi_ini_s
-                     zsm_i_ini(ji,jj) = rn_smi_ini_s
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_s
+                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
+                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
+                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
+                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
+                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
                   ENDIF
                END DO
@@ -181,5 +187,4 @@
          zh_i_ini(:,:,:) = 0._wp 
          za_i_ini(:,:,:) = 0._wp
-         zv_i_ini(:,:,:) = 0._wp
 
          DO jj = 1, jpj
@@ -188,128 +193,109 @@
                IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
 
-                  ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
-!                  ztests  = 0 
-
-                  DO i_fill = jpl, 1, -1
-
-!                     IF( ztests .NE. 4 ) THEN
-                     IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
-                        !----------------------------
-                        ! fill the i_fill categories
-                        !----------------------------
-                        ! *** 1 category to fill
-                        IF ( i_fill .EQ. 1 ) THEN
-                           zh_i_ini(ji,jj,    1)   = zht_i_ini(ji,jj)
-                           za_i_ini(ji,jj,    1)   = zat_i_ini(ji,jj)
-                           zh_i_ini(ji,jj,2:jpl)   = 0._wp
-                           za_i_ini(ji,jj,2:jpl)   = 0._wp
-                        ELSE
-
-                           ! *** >1 categores to fill
-                           !--- Ice thicknesses in the i_fill - 1 first categories
-                           DO jl = 1, i_fill - 1
-                              zh_i_ini(ji,jj,jl) = hi_mean(jl)
-                           END DO
+                  !--- jl0: most likely index where cc will be maximum
+                  jl0 = jpl
+                  DO jl = 1, jpl
+                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
+                        jl0 = jl
+                        CYCLE
+                     ENDIF
+                  END DO
+
+                  ! 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
+                     zh_i_ini(ji,jj,:) = 0._wp 
+                     za_i_ini(ji,jj,:) = 0._wp
+                     itest(:) = 0
+
+                     ! *** case very thin ice: fill only category 1
+                     IF ( i_fill == 1 ) THEN
+                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
+                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
+
+                     ! *** case ice is thicker: fill categories >1
+                     ELSE
+
+                        ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
+                        DO jl = 1, i_fill-1
+                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
+                        END DO
                
-                           !--- jl0: most likely index where cc will be maximum
-                           DO jl = 1, jpl
-                              IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. &
-                                 & ( zht_i_ini(ji,jj) <= hi_max(jl)   ) ) THEN
-                                 jl0 = jl
-                              ENDIF
-                           END DO
-                           jl0 = MIN(jl0, i_fill)
-               
-                           !--- Concentrations
-                           za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
-                           DO jl = 1, i_fill - 1
-                              IF( jl .NE. jl0 )THEN
-                                 zsigma             = 0.5 * zht_i_ini(ji,jj)
-                                 zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma
-                                 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
-                              ENDIF
-                           END DO
-               
-                           zA = 0. ! sum of the areas in the jpl categories 
-                           DO jl = 1, i_fill - 1
-                              zA = zA + za_i_ini(ji,jj,jl)
-                           END DO
-                           za_i_ini(ji,jj,i_fill)   = zat_i_ini(ji,jj) - zA ! ice conc in the last category
-                           IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
-         
-                           !--- Ice thickness in the last category
-                           zV = 0. ! sum of the volumes of the N-1 categories
-                           DO jl = 1, i_fill - 1
-                              zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl)
-                           END DO
-                           zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 
-                           IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
-
-                           !--- volumes
-                           zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:)
-                           IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
-
-                        ENDIF ! i_fill
-
-                        !---------------------
-                        ! Compatibility tests
-                        !---------------------
-                        ! Test 1: area conservation
-                        zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons )
-                        IF ( zconv .LT. 1.0e-6 ) THEN
-                           ztest_1 = 1
-                        ELSE 
-                          ztest_1 = 0
-                        ENDIF
-
-                        ! Test 2: volume conservation
-                        zV_cons = SUM(zv_i_ini(ji,jj,:))
-                        zconv = ABS(zvt_i_ini(ji,jj) - zV_cons)
-
-                        IF( zconv .LT. 1.0e-6 ) THEN
-                           ztest_2 = 1
-                        ELSE
-                           ztest_2 = 0
-                        ENDIF
-
-                        ! Test 3: thickness of the last category is in-bounds ?
-                        IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN
-                           ztest_3 = 1
-                        ELSE
-                           ztest_3 = 0
-                        ENDIF
-
-                        ! Test 4: positivity of ice concentrations
-                        ztest_4 = 1
-                        DO jl = 1, jpl
-                           IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 
-                              ztest_4 = 0
+                        !--- Concentrations
+                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
+                        DO jl = 1, i_fill - 1
+                           IF( jl /= jl0 )THEN
+                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
+                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
                            ENDIF
                         END DO
-
-                     ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
- 
-                     ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
-
-                  END DO ! i_fill
-
-                  IF(lwp) THEN 
-                     WRITE(numout,*) ' ztests : ', ztests
-                     IF( ztests .NE. 4 )THEN
-                        WRITE(numout,*)
-                        WRITE(numout,*) ' !!!! ALERT                  !!! '
-                        WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
-                        WRITE(numout,*)
-                        WRITE(numout,*) ' *** ztests is not equal to 4 '
-                        WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
-                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
-                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
-                     ENDIF ! ztests .NE. 4
+               
+                        ! Concentration in the last (i_fill) category
+                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
+
+                        ! Ice thickness in the last (i_fill) category
+                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
+                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
+
+                        ! clem: correction if concentration of upper cat is greater than lower cat
+                        !       (it should be a gaussian around jl0 but sometimes it is not)
+                        IF ( jl0 /= jpl ) THEN
+                           DO jl = jpl, jl0+1, -1
+                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
+                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
+                                 zh_i_ini(ji,jj,jl    ) = 0._wp
+                                 za_i_ini(ji,jj,jl    ) = 0._wp
+                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
+                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
+                              END IF
+                           ENDDO
+                        ENDIF
+
+                     ENDIF ! case ice is thick or thin
+
+                     !---------------------
+                     ! Compatibility tests
+                     !---------------------
+                     ! Test 1: area conservation
+                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
+                     IF ( zconv < epsi06 ) itest(1) = 1
+                     
+                     ! Test 2: volume conservation
+                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &
+                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
+                     IF ( zconv < epsi06 ) itest(2) = 1
+                     
+                     ! Test 3: thickness of the last category is in-bounds ?
+                     IF ( zh_i_ini(ji,jj,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_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
+                     END DO
+                     !                                      !============================
+                  END DO                                    ! end iteration on categories
+                  !                                         !============================
+
+                  IF( lwp .AND. SUM(itest) /= 4 ) THEN 
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
+                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
+                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
+                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
                   ENDIF
-      
-               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp
-
-            ENDDO   
-         ENDDO   
+               
+               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
+            
+            ENDDO
+         ENDDO
 
          !---------------------------------------------------------------------
@@ -468,8 +454,13 @@
       sxyage (:,:,:)  = 0._wp
 
-
-      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini )
+!!!clem
+!!      ! Output the initial state and forcings
+!!      CALL dia_wri_state( 'output.init', nit000 )
+!!!
+      
+      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini )
       CALL wrk_dealloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
       CALL wrk_dealloc( jpi, jpj,      zswitch )
+      Call wrk_dealloc( 4,             itest )
 
    END SUBROUTINE lim_istate
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6801)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6853)
@@ -38,4 +38,6 @@
    USE bdyice_lim
 #endif
+   !!!clem
+!!   USE iom
 
    IMPLICIT NONE
@@ -140,9 +142,10 @@
                                                                              !   ocean surface (ssh_m) if ice is not embedded
                                                                              !   ice top surface if ice is embedded   
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitch1, zswitch2              ! dummy arrays
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
       REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice
 
       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
+      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
       REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
       !!-------------------------------------------------------------------
@@ -152,5 +155,5 @@
       CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
       CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
-      CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf )
+      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
 
 #if defined key_agrif 
@@ -294,7 +297,11 @@
             zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
 
+            ! masks
+            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
             ! switches
-            zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )
-            zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )
+            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
+            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
 
          END DO
@@ -436,13 +443,14 @@
                   
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
-                     &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &           ) * zswitch2(ji,jj)
+                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
                   ! Bouillon 2013
                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
                   !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
+                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
 
                END DO
@@ -480,13 +488,14 @@
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
-                     &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &           ) * zswitch1(ji,jj)
+                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
+                     &           ) * zmaskU(ji,jj)
                   ! Bouillon 2013
                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
                   !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
+                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
                END DO
             END DO
@@ -525,13 +534,14 @@
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
-                     &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast
-                     &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &           ) * zswitch1(ji,jj)
+                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskU(ji,jj)
                   ! Bouillon 2013
                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
                   !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
+                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
                END DO
             END DO
@@ -564,17 +574,18 @@
                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
 
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
-                     &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast
-                     &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &           ) * zswitch2(ji,jj)
+                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
                   ! Bouillon 2013
                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
                   !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
+                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
                END DO
             END DO
@@ -603,29 +614,5 @@
       !
       !------------------------------------------------------------------------------!
-      ! 4) Limit ice velocities when ice is thin
-      !------------------------------------------------------------------------------!
-      ! 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 = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            rswitch      = MAX( 0._wp, SIGN( 1._wp, vt_i(ji,jj) - zvmin ) )
-            u_ice(ji,jj) = ( rswitch * u_ice(ji,jj) + (1._wp - rswitch) * u_oce(ji,jj) ) * zswitch1(ji,jj)
-            v_ice(ji,jj) = ( rswitch * v_ice(ji,jj) + (1._wp - rswitch) * v_oce(ji,jj) ) * zswitch2(ji,jj)
-         END DO
-      END DO
-
-      CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )
-
-#if defined key_agrif
-      CALL agrif_interp_lim3( 'U' )
-      CALL agrif_interp_lim3( 'V' )
-#endif
-#if defined key_bdy
-      CALL bdy_ice_lim_dyn( 'U' )
-      CALL bdy_ice_lim_dyn( 'V' )
-#endif         
-
-      !------------------------------------------------------------------------------!
-      ! 5) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
       !------------------------------------------------------------------------------!
       DO jj = 1, jpjm1
@@ -665,5 +652,6 @@
             ! delta at T points
             zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )  
-            delta_i(ji,jj) = zdelta + rn_creepl
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
 
          END DO
@@ -678,5 +666,5 @@
 
       !------------------------------------------------------------------------------!
-      ! 6) Control prints of residual and charge ellipse
+      ! 5) Control prints of residual and charge ellipse
       !------------------------------------------------------------------------------!
       !
@@ -713,9 +701,19 @@
       ENDIF
       !
+
+!!      ! velocity
+!!      IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN 
+!!         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component
+!!         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component
+!!      ENDIF
+!!      IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i    )        ! area for categories
+!!      IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i   )        ! thickness for categories
+     
+      
       CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
       CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
       CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
       CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
-      CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf )
+      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
 
    END SUBROUTINE lim_rhg
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6801)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6853)
@@ -636,10 +636,10 @@
       INTEGER  :: ji, jk, jl             ! dummy loop indices
       INTEGER  :: ijpij, i_fill, jl0  
-      REAL(wp) :: zarg, zV, zconv, zdh
+      REAL(wp) :: zarg, zV, zconv, zdh, zdv
       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 )
+      Call wrk_alloc( 4, itest )
       !--------------------------------------------------------------------
       ! initialisation of variables
@@ -657,80 +657,95 @@
          IF( zhti(ji) > 0._wp ) THEN
 
-         ! initialisation of tests
-         itest(:)  = 0
+            ! find which category (jl0) the input ice thickness falls into
+            jl0 = jpl
+            DO jl = 1, jpl
+               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
+                  jl0 = jl
+                  CYCLE
+               ENDIF
+            END DO
+
+            ! 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
+            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
+               itest(:)        = 0           
+               
+               ! *** 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 in the (i_fill-1) cat by hmean 
+                  DO jl = 1, i_fill - 1
+                     zht_i(ji,jl) = hi_mean(jl)
+                  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 ) THEN
+                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
+                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
+                     ENDIF
+                  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 ) / MAX( za_i(ji,i_fill), epsi10 ) 
+                  
+                  ! clem: correction if concentration of upper cat is greater than lower cat
+                  !       (it should be a gaussian around jl0 but sometimes it is not)
+                  IF ( jl0 /= jpl ) THEN
+                     DO jl = jpl, jl0+1, -1
+                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
+                           zdv = zht_i(ji,jl) * za_i(ji,jl)
+                           zht_i(ji,jl    ) = 0._wp
+                           za_i (ji,jl    ) = 0._wp
+                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
+                        END IF
+                     ENDDO
+                  ENDIF
+               
+               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
             
-            ! 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
+               ! 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
                
-               ! find which category (jl0) the input ice thickness falls into
-               jl0 = i_fill
+               ! 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 ( ( 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
-                                                           !============================
+                  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
Index: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6801)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6853)
@@ -96,10 +96,10 @@
          END DO
          CALL lbc_lnk( z2d, 'T', 1. )
-         CALL iom_put( "uice_ipa"     , u_ice * zswi     )       ! ice velocity u component
-         CALL iom_put( "vice_ipa"     , v_ice * zswi     )       ! ice velocity v component
-         CALL iom_put( "icevel"       , z2d   * zswi     )       ! ice velocity module
+         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component
+         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component
+         CALL iom_put( "icevel"       , z2d        )       ! ice velocity module
       ENDIF
 
-      IF ( iom_use( "tau_icebfr" ) )    CALL iom_put( "tau_icebfr"  , tau_icebfr * zswi      )  ! ice friction with ocean bottom (landfast ice)  
+      IF ( iom_use( "tau_icebfr" ) )    CALL iom_put( "tau_icebfr"  , tau_icebfr             )  ! ice friction with ocean bottom (landfast ice)  
       !
       IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age
@@ -117,12 +117,12 @@
       CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content
       CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume
-      CALL iom_put( "utau_ice"    , utau_ice  * zswi    )        ! wind stress over ice along i-axis at I-point
-      CALL iom_put( "vtau_ice"    , vtau_ice  * zswi    )        ! wind stress over ice along j-axis at I-point
+      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point
+      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point
       CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation 
       CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
 
       CALL iom_put( "icestr"      , strength * zswi )    ! ice strength
-      CALL iom_put( "idive"       , divu_i * 1.0e8   * zswi )    ! divergence
-      CALL iom_put( "ishear"      , shear_i * 1.0e8  * zswi )    ! shear
+      CALL iom_put( "idive"       , divu_i * 1.0e8      )    ! divergence
+      CALL iom_put( "ishear"      , shear_i * 1.0e8     )    ! shear
       CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume
       
@@ -232,78 +232,100 @@
       !!   4.0  !  2013-06  (C. Rousset)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt               ! ocean time-step index)
-      INTEGER, INTENT( in ) ::   kid , kh_i       
+      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index)
+      INTEGER, INTENT( in )   ::   kid , kh_i
+      INTEGER                 ::   nz_i, jl
+      REAL(wp), DIMENSION(jpl) :: jcat
       !!----------------------------------------------------------------------
-
-      CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2" ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" ,   &
+      DO jl = 1, jpl
+         jcat(jl) = REAL(jl)
+      ENDDO
+      
+      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up")
+
+      CALL histdef( kid, "sithic", "Ice thickness"           , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "siconc", "Ice concentration"       , "%"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sitemp", "Ice temperature"         , "C"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sivelu", "i-Ice speed "            , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sivelv", "j-Ice speed "            , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sisflx", "Solar flux over ocean"     , "w/m2" ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" ,   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    ,   &
+      CALL histdef( kid, "sisali", "Ice salinity"            , "PSU"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sivolu", "Ice volume"              , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+
+      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxdyn", "Ice dynamic production"  , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxopw", "Ice open water prod"     , "m/s"    ,   &
       &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""       ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxsni", "Snow ice production "    , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxres", "Ice prod from limupdate" , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxbom", "Ice bottom melt"         , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxsum", "Ice surface melt"        , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+
+      CALL histdef( kid, "sithicat", "Ice thickness"         , "m"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "siconcat", "Ice concentration"     , "%"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sisalcat", "Ice salinity"           , ""      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sitemcat", "Ice temperature"       , "C"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "snthicat", "Snw thickness"         , "m"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sntemcat", "Snw temperature"       , "C"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
 
       CALL histend( kid, snc4set )   ! end of the file definition
 
-      CALL histwrite( kid, "iicethic", kt, htm_i         , jpi*jpj, (/1/) )    
-      CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sithic", kt, htm_i         , jpi*jpj, (/1/) )    
+      CALL histwrite( kid, "siconc", kt, at_i          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sitemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivelu", kt, u_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivelv", kt, v_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sistru", kt, utau_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sistrv", kt, vtau_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) )
       CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesali", kt, smt_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevolu", kt, vt_i           , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
-
-      CALL histwrite( kid, "iicebopr", kt, wfx_bog        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedypr", kt, wfx_dyn        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicelapr", kt, wfx_opw        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesipr", kt, wfx_sni        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicerepr", kt, wfx_res        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sisali", kt, smt_i          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivolu", kt, vt_i           , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
+
+      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxdyn", kt, wfx_dyn        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxopw", kt, wfx_opw        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxsni", kt, wfx_sni        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxres", kt, wfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) )
+
+      CALL histwrite( kid, "sithicat", kt, ht_i        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "siconcat", kt, a_i         , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sisalcat", kt, sm_i        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sitemcat", kt, tm_i - rt0  , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "snthicat", kt, ht_s        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )    
 
       ! Close the file
