Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4332)
@@ -292,6 +292,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dh_i_surf2D, dh_i_bott2D, fstbif, fsup2D, focea2D, q_s
 
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   patho_case ! number of the pathological case (if any, of course)
-
    !!--------------------------------------------------------------------------
    !! * Ice global state variables
@@ -396,6 +394,6 @@
    LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F)
    LOGICAL               , PUBLIC ::   ln_nicep      = .FALSE.             !: flag for sea-ice points output (T) or not (F)
-   REAL(wp)              , PUBLIC ::   cai           = 1.40e-3            !: atmospheric drag over sea ice
-   REAL(wp)              , PUBLIC ::   cao           = 1.00e-3            !: atmospheric drag over ocean
+   REAL(wp)              , PUBLIC ::   cai           = 1.4e-3            !: atmospheric drag over sea ice
+   REAL(wp)              , PUBLIC ::   cao           = 1.0e-3            !: atmospheric drag over ocean
    REAL(wp)              , PUBLIC ::   amax          = 0.99               !: maximum ice concentration
    !
@@ -465,7 +463,4 @@
          &      fsup2D     (jpi,jpj) , focea2D    (jpi,jpj) , q_s   (jpi,jpj) , STAT=ierr(ii) )
 
-      ii = ii + 1
-      ALLOCATE( patho_case(jpi, jpj, jpl)  , STAT=ierr(ii) )
-
       ! * Ice global state variables
       ii = ii + 1
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4332)
@@ -88,5 +88,5 @@
       !
       !                                ! Initial sea-ice state
-      IF( .NOT.ln_rstart ) THEN              ! start from rest
+      IF( .NOT. ln_rstart ) THEN              ! start from rest
          numit = 0
          numit = nit000 - 1
@@ -101,8 +101,10 @@
       ENDIF
       !
-      CALL lim_sbc_init                ! ice surface boundary condition   
-      !
-      fr_i(:,:) = at_i(:,:)           ! initialisation of sea-ice fraction
-      tn_ice(:,:,:) = t_su(:,:,:)
+      hi_max(jpl) = 99._wp  ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl)
+      !
+      CALL lim_sbc_init                 ! ice surface boundary condition   
+      !
+      fr_i(:,:) = at_i(:,:)             ! initialisation of sea-ice fraction
+      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
       !
       nstart = numit  + nn_fsbc      
@@ -132,8 +134,8 @@
       READ  ( numnam_ice , namicerun )
       !
-      IF( lk_mpp .AND. ln_nicep ) THEN
-         ln_nicep = .FALSE.
-         CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' )
-      ENDIF
+      !IF( lk_mpp .AND. ln_nicep ) THEN
+      !   ln_nicep = .FALSE.
+      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' )
+      !ENDIF
       !
       IF(lwp) THEN                        ! control print
@@ -159,4 +161,6 @@
       !! ** Purpose :   Initializes the ice thickness distribution
       !! ** Method  :   ...
+      !!    Note    : hi_max(jpl) is here set up to a value close to 7 m for
+      !!              limistate (only) and is changed to 99 m in ice_init
       !!------------------------------------------------------------------
       INTEGER  ::   jl, jm               ! dummy loop index
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90	(revision 4332)
@@ -44,6 +44,8 @@
       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
-      REAL(wp) :: epsi06 = 1.0e-6
-  
+      REAL(wp)                     ::   epsi06 = 1.0e-6
+      REAL(wp)                     ::   zc1, zc2, zc3, zx1, zdh   ! local scalars
+      REAL(wp), DIMENSION(0:jpl)   ::   zhi_max         !:Boundary of ice thickness categories in thickness space
+ 
       IF( nn_timing == 1 )  CALL timing_start('limcat_1D')
       !--------------------------------------------------------------------
@@ -51,7 +53,7 @@
       !--------------------------------------------------------------------
       ijpij = SIZE(zhti,1)
-      zht_i(1:ijpij,1:jpl) = 0.d0
-      zht_s(1:ijpij,1:jpl) = 0.d0
-      za_i (1:ijpij,1:jpl) = 0.d0
+      zht_i(1:ijpij,1:jpl) = 0._wp
+      zht_s(1:ijpij,1:jpl) = 0._wp
+      za_i (1:ijpij,1:jpl) = 0._wp
 
       !------------------------------------------------------------------------------------
@@ -66,4 +68,16 @@
       !         fulfills ice volume concervation between input and output (ztests=4) 
       !--------------------------------------------------------------------------------------
+
+      !- Thickness categories boundaries
+      !    hi_max is calculated in iceini.F90 but since limcat_1D.F90 routine 
+      !    is called before (in bdydta.F90), one must recalculate it.   
+      !    Note clem: there may be a way of doing things cleaner
+      !----------------------------------
+      zhi_max(:) = 0._wp
+      zc1 =  3._wp / REAL( jpl , wp ) ; zc2 = 10._wp * zc1 ; zc3 =  3._wp
+      DO jl = 1, jpl
+         zx1 = REAL( jl-1 , wp ) / REAL( jpl , wp )
+         zhi_max(jl) = zhi_max(jl-1) + zc1 + zc2 * ( 1._wp + TANH( zc3 * ( zx1 - 1._wp ) ) )
+      END DO
  
       ! ----------------------------------------
@@ -71,7 +85,7 @@
       ! ----------------------------------------
       DO ji = 1, ijpij
-        ! snow thickness in each category 
-         zht_s(ji,1:jpl) = zhts(ji)
          
+         IF( zhti(ji) > 0._wp ) THEN
+
          ! initialisation of tests
          ztest_1 = 0
@@ -87,6 +101,6 @@
             
             ! initialisation of ice variables for each try
-            zht_i(ji,1:jpl) = 0.d0
-            za_i (ji,1:jpl) = 0.d0
+            zht_i(ji,1:jpl) = 0._wp
+            za_i (ji,1:jpl) = 0._wp
             
             ! *** case very thin ice: fill only category 1
@@ -94,5 +108,4 @@
                zht_i(ji,1) = zhti(ji)
                za_i (ji,1) = zai (ji)
-               
                ! *** case ice is thicker: fill categories >1
             ELSE
@@ -100,5 +113,5 @@
                ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2 
                DO jl = 1, i_fill - 1
-                  zht_i(ji,jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.
+                  zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) * 0.5_wp
                END DO
                
@@ -106,5 +119,5 @@
                jl0 = i_fill
                DO jl = 1, i_fill
-                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
+                  IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_max(jl) ) ) THEN
                      jl0 = jl
 		     CYCLE
@@ -116,5 +129,5 @@
                DO jl = 1, i_fill - 1
                   IF ( jl == jl0 ) CYCLE
-                  zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) / 2. )
+                  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
@@ -141,10 +154,10 @@
             
             ! Test 3: thickness of the last category is in-bounds ?
-            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1
+            IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1
             
             ! Test 4: positivity of ice concentrations
             ztest_4 = 1
             DO jl = 1, i_fill
-               IF ( za_i(ji,jl) < 0.0d0 ) ztest_4 = 0
+               IF ( za_i(ji,jl) < 0._wp ) ztest_4 = 0
             END DO
             
@@ -164,6 +177,25 @@
          !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
          !ENDIF
-            
+         
+         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( zhi_max(jl), zht_i(ji,jl) + zdh )
+               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
+            ENDIF
+         ENDDO
+      ENDDO
 
       IF( nn_timing == 1 )  CALL timing_stop('limcat_1D')
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90	(revision 4220)
+++ 	(revision )
@@ -1,180 +1,0 @@
-MODULE limcat_2D
-   !!======================================================================
-   !!                     ***  MODULE  limcat_2D  ***
-   !!  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      :  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
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! Accessibility
-   PUBLIC lim_cat_2D      
-
-CONTAINS
-
-   SUBROUTINE lim_cat_2D(zhti,zhts,zai,zht_i,zht_s,za_i)
-      !! Local variables
-      INTEGER  :: ji, jj, jk, jl             ! dummy loop indices
-      INTEGER  :: ijpi, ijpj, i_fill, jl0, ztest_1, ztest_2, ztest_3, ztest_4, ztests  
-      REAL(wp) :: zarg, zV, zconv
-      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
-      REAL(wp) :: epsi06 = 1.0e-6
-  
-      !--------------------------------------------------------------------
-      ! initialisation of variables
-      !--------------------------------------------------------------------
-      ijpi = SIZE(zhti,1)
-      ijpj = SIZE(zhti,2)
-      zht_i(1:ijpi,1:ijpj,1:jpl) = 0.d0
-      zht_s(1:ijpi,1:ijpj,1:jpl) = 0.d0
-      za_i (1:ijpi,1:ijpj,1:jpl) = 0.d0
-
-      !------------------------------------------------------------------------------------
-      ! 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 (ztests=4) 
-      !--------------------------------------------------------------------------------------
- 
-      ! ----------------------------------------
-      ! distribution over the jpl ice categories
-      ! ----------------------------------------
-      DO jj = 1, ijpj
-         DO ji = 1, ijpi
-
-            ! snow thickness in each category 
-            zht_s(ji,jj,1:jpl) = zhts(ji,jj)
-
-            ! initialisation of tests
-            ztest_1 = 0
-            ztest_2 = 0
-            ztest_3 = 0
-            ztest_4 = 0
-            ztests  = 0
-         
-            i_fill = jpl + 1                                    !====================================
-            DO WHILE ( ( ztests /= 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,jj,1:jpl) = 0.d0
-               za_i (ji,jj,1:jpl) = 0.d0
-
-               ! *** case very thin ice: fill only category 1
-               IF ( i_fill == 1 ) THEN
-                  zht_i(ji,jj,1) = zhti(ji,jj)
-                  za_i (ji,jj,1) = zai (ji,jj)
-            
-               ! *** case ice is thicker: fill categories >1
-               ELSE
-
-                  ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2 
-                  DO jl = 1, i_fill - 1
-                     zht_i(ji,jj,jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.
-                  END DO
-
-                  ! find which category (jl0) the input ice thickness falls into
-                  jl0 = i_fill
-                  DO jl = 1, i_fill
-                     IF ( ( zhti(ji,jj) >= hi_max(jl-1) ) .AND. ( zhti(ji,jj) < hi_max(jl) ) ) THEN
-                        jl0 = jl
-                        CYCLE
-                     ENDIF
-                  END DO
-
-                  ! Concentrations in the (i_fill-1) categories 
-                  za_i(ji,jj,jl0) = zai(ji,jj) / SQRT(REAL(jpl))
-                  DO jl = 1, i_fill - 1
-                     IF ( jl == jl0 ) CYCLE
-                     zarg           = ( zht_i(ji,jj,jl) - zhti(ji,jj) ) / ( zhti(ji,jj) / 2. )
-                     za_i(ji,jj,jl) =   za_i (ji,jj,jl0) * EXP(-zarg**2)
-                  END DO 
-
-                  ! Concentration in the last (i_fill) category
-                  za_i(ji,jj,i_fill) = zai(ji,jj) - SUM( za_i(ji,jj,1:i_fill-1) )
-
-                  ! Ice thickness in the last (i_fill) category
-                  zV = SUM( za_i(ji,jj,1:i_fill-1) * zht_i(ji,jj,1:i_fill-1) )
-                  zht_i(ji,jj,i_fill) = ( zhti(ji,jj)*zai(ji,jj) -  zV ) / za_i(ji,jj,i_fill) 
- 
-               ENDIF ! case ice is thick or thin
-
-               !---------------------
-               ! Compatibility tests
-               !--------------------- 
-               ! Test 1: area conservation
-               zconv = ABS( zai(ji,jj) - SUM( za_i(ji,jj,1:jpl) ) )
-               IF ( zconv < epsi06 ) ztest_1 = 1
- 
-               ! Test 2: volume conservation
-               zconv = ABS( zhti(ji,jj)*zai(ji,jj) - SUM( za_i(ji,jj,1:jpl)*zht_i(ji,jj,1:jpl) ) )
-               IF ( zconv < epsi06 ) ztest_2 = 1
- 
-               ! Test 3: thickness of the last category is in-bounds ?
-               IF ( zht_i(ji,jj,i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1
- 
-               ! Test 4: positivity of ice concentrations
-               ztest_4 = 1
-               DO jl = 1, i_fill
-                  IF ( za_i(ji,jj,jl) < 0.0d0 ) ztest_4 = 0
-               END DO
-   
-               ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
-                                                              !============================
-            END DO                                            ! end iteration on categories
-                                                              !============================
-            ! Check if tests have passed (i.e. volume conservation...)
-            !IF ( ztests .NE. 4 ) THEN
-            !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
-            !   WRITE(numout,*) ' !! ALERT categories distribution !!'
-            !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
-            !   WRITE(numout,*) ' *** ztests is not equal to 4 '
-            !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
-            !ENDIF
-
-         END DO ! i loop
-      END DO ! j loop
-
-   END SUBROUTINE lim_cat_2D
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :         Empty module          NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_cat_2D          ! Empty routine
-   END SUBROUTINE lim_cat_2D
-#endif
-
-   !!======================================================================
-END MODULE limcat_2D
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 4332)
@@ -29,5 +29,5 @@
 
    LOGICAL  ::   linit = .TRUE.              ! initialization flag (set to flase after the 1st call)
-   REAL(wp) ::   epsi04 = 1e-04              ! constant
+   REAL(wp) ::   epsi04 = 1.e-04              ! constant
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient
 
@@ -54,7 +54,7 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied
       !
-      INTEGER  ::  ji, jj            ! dummy loop indices
-      INTEGER  ::  its, iter, ierr   ! local integers
-      REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! local scalars
+      INTEGER  ::  ji, jj                   ! dummy loop indices
+      INTEGER  ::  its, iter, ierr          ! local integers
+      REAL(wp) ::   zalfa, zrlxint, zconv   ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:) ::   zrlx, zflu, zflv, zdiv0, zdiv, ztab0
       CHARACTER(lc) ::   charout   ! local character
@@ -79,5 +79,4 @@
       zalfa = 0.5_wp                      ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
       its   = 100                         ! Maximum number of iteration
-      zeps  =  2._wp * epsi04
       !
       ztab0(:, : ) = ptab(:,:)      ! Arrays initialization
@@ -94,5 +93,5 @@
       iter  = 0
       !
-      DO WHILE( zconv > zeps .AND. iter <= its )   ! Sub-time step loop
+      DO WHILE( zconv > ( 2._wp * epsi04 ) .AND. iter <= its )   ! Sub-time step loop
          !
          iter = iter + 1                                 ! incrementation of the sub-time step number
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4332)
@@ -76,4 +76,7 @@
       !!                6) Lateral boundary conditions
       !!
+      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
+      !!              where there is no ice (clem: I do not know why but it is mandatory) 
+      !!
       !! History :
       !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
@@ -84,5 +87,5 @@
       !! * Local variables
       INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
-      REAL(wp)   :: epsi06, epsi20, ztmelts
+      REAL(wp)   :: epsi20, 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
@@ -98,5 +101,4 @@
       CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini )
 
-      epsi06   = 1.0e-6
       epsi20   = 1.0e-20
       IF(lwp) WRITE(numout,*)
@@ -308,6 +310,5 @@
       !---------------------------------------------------------------------
 
-      ! Ice concentration, thickness and volume, snow depth, ice
-      ! salinity, ice age, surface temperature
+      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
       DO jl = 1, jpl ! loop over categories
          DO jj = 1, jpj
@@ -315,10 +316,18 @@
                a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
                ht_i(ji,jj,jl)  = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj))  ! ice thickness
-               ht_s(ji,jj,jl)  = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj))     ! snow depth
-               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
+               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zhm_s_ini( zhemis(ji,jj) ) / zhm_i_ini( zhemis(ji,jj) ) )  ! snow depth
+               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ! salinity
                o_i(ji,jj,jl)   = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age
-               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp
- 
-               ! ice volume, snow volume, salt content, age content
+               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * 270.0 ! surf temp
+
+               ! This case below should not be used if (ht_s/ht_i) is ok in namelist
+               ! 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 * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
+               ! recompute ht_i, ht_s avoiding out of bounds values
+               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
+               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn )
+
+               ! ice volume, salt content, age content
                v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
                v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4332)
@@ -43,5 +43,5 @@
    PUBLIC   lim_itd_me_alloc        ! called by iceini.F90
 
-   REAL(wp) ::   epsi11 = 1.e-11_wp   ! constant values
+   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values
    REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values
    REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values
@@ -50,20 +50,18 @@
    ! Variables shared among ridging subroutines
    !-----------------------------------------------------------------------
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
-   !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/
-   !                                                           !  closing associated w/ category n
-   !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/
+   !                                                     ! closing associated w/ category n
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
 
    REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
    REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer
-   REAL(wp), PARAMETER ::   kamax   = 1.0
+   REAL(wp), PARAMETER ::   kamax   = 1.0_wp    ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
 
    REAL(wp) ::   Cp                             ! 
@@ -183,8 +181,4 @@
       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
       !-----------------------------------------------------------------------------!
-      ! Set hi_max(ncat) to a big value to ensure that all ridged ice is thinner than hi_max(ncat).
-
-      hi_max(jpl) = 999.99
-
       Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE
       !
@@ -265,5 +259,5 @@
                ! Reduce the closing rate if more than 100% of the open water 
                ! would be removed.  Reduce the opening rate proportionately.
-               IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN
+               IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN
                   w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
                   IF ( w1 .GT. ato_i(ji,jj)) THEN
@@ -285,5 +279,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF ( a_i(ji,jj,jl) > epsi11 .AND. athorn(ji,jj,jl) > 0._wp )THEN
+                  IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN
                      w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
                      IF ( w1  >  a_i(ji,jj,jl) ) THEN
@@ -316,5 +310,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN
+               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN
                   closing_net(ji,jj) = 0._wp
                   opning     (ji,jj) = 0._wp
@@ -358,5 +352,5 @@
          DO ji = 1, jpi
 
-            IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true.
+            IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true.
 
             dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
@@ -377,5 +371,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( ABS( asum(ji,jj) - kamax)  >  epsi11 ) THEN   ! there is a bug
+            IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
                WRITE(numout,*) ' '
                WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
@@ -400,7 +394,6 @@
 
       !-----------------------------------------------------------------------------!
-      ! 6) Updating state variables and trend terms
+      ! 6) Updating state variables and trend terms (done in limupdate)
       !-----------------------------------------------------------------------------!
-
       CALL lim_var_glo2eqv
       CALL lim_itd_me_zapsmall
@@ -419,19 +412,4 @@
          END DO
       END DO
-
-      !-----------------
-      !  Trend terms
-      !-----------------
-
-      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:)
-      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:)
-      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:)
-      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)  
-      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)   
-      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)  
-      d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:)
-      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:)
-      d_smv_i_trp(:,:,:)   = 0._wp
-      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:)
 
       IF(ln_ctl) THEN     ! Control print
@@ -491,64 +469,4 @@
       ! -------------------------------
 
-      !-------------------------!
-      ! Back to initial values
-      !-------------------------!
-
-      ! update of fields will be made later in lim update
-      u_ice(:,:)   = old_u_ice(:,:)
-      v_ice(:,:)   = old_v_ice(:,:)
-      a_i(:,:,:)   = old_a_i(:,:,:)
-      v_s(:,:,:)   = old_v_s(:,:,:)
-      v_i(:,:,:)   = old_v_i(:,:,:)
-      e_s(:,:,:,:) = old_e_s(:,:,:,:)
-      e_i(:,:,:,:) = old_e_i(:,:,:,:)
-      oa_i(:,:,:)  = old_oa_i(:,:,:)
-      IF(  num_sal == 2  )   smv_i(:,:,:) = old_smv_i(:,:,:)
-
-      !----------------------------------------------------!
-      ! Advection of ice in a free cell, newly ridged ice
-      !----------------------------------------------------!
-
-      ! to allow for thermodynamics to melt new ice
-      ! we immediately advect ice in free cells
-
-      ! heat content has to be corrected before ice volume
-!clem@order
-!      DO jl = 1, jpl
-!         DO jk = 1, nlay_i
-!            DO jj = 1, jpj
-!               DO ji = 1, jpi
-!                  IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. &
-!                     ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN
-!                     old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl)
-!                     d_e_i_trp(ji,jj,jk,jl) = 0._wp
-!                  ENDIF
-!               END DO
-!            END DO
-!         END DO
-!      END DO
-!
-!      DO jl = 1, jpl
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. &
-!                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN
-!                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl)
-!                  d_v_i_trp (ji,jj,jl)   = 0._wp
-!                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl)
-!                  d_a_i_trp (ji,jj,jl)   = 0._wp
-!                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl)
-!                  d_v_s_trp (ji,jj,jl)   = 0._wp
-!                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl)
-!                  d_e_s_trp (ji,jj,1,jl) = 0._wp
-!                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl)
-!                  d_oa_i_trp(ji,jj,jl)   = 0._wp
-!                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl)
-!                  d_smv_i_trp(ji,jj,jl)  = 0._wp
-!               ENDIF
-!            END DO
-!         END DO
-!      END DO
-!clem@order
       ENDIF  ! ln_limdyn=.true.
       !
@@ -605,6 +523,5 @@
                DO ji = 1, jpi
                   !
-                  IF(  a_i(ji,jj,jl)    > epsi11  .AND.     &
-                       athorn(ji,jj,jl) > 0._wp   ) THEN
+                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
                      hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
                      !----------------------------
@@ -624,5 +541,5 @@
                         * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )   
 !!gm Optimization:  (a**3-b**3)/(a-b) = a*a+ab+b*b   ==> less costly operations even if a**3 is replaced by a*a*a...                    
-                  ENDIF            ! aicen > epsi11
+                  ENDIF            ! aicen > epsi10
                   !
                END DO ! ji
@@ -681,5 +598,5 @@
          DO jj = 2, jpj - 1
             DO ji = 2, jpi - 1
-               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is
+               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is
                   ! present
                   zworka(ji,jj) = 4.0 * strength(ji,jj)              &
@@ -721,5 +638,5 @@
          DO jj = 1, jpj - 1
             DO ji = 1, jpi - 1
-               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN       ! ice is present
+               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present
                   numts_rm = 1 ! number of time steps for the running mean
                   IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
@@ -794,5 +711,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( ato_i(ji,jj) > epsi11 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj)
+            IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj)
             ELSE                               ;   Gsum(ji,jj,0) = 0._wp
             ENDIF
@@ -804,5 +721,5 @@
          DO jj = 1, jpj 
             DO ji = 1, jpi
-               IF( a_i(ji,jj,jl) .GT. epsi11 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
+               IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
                ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
                ENDIF
@@ -887,10 +804,9 @@
       IF ( raftswi == 1 ) THEN
 
-         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN
+         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN
             DO jl = 1, jpl
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. &
-                        epsi11 ) THEN
+                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN
                         WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
                         WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
@@ -938,5 +854,5 @@
             DO ji = 1, jpi
 
-               IF (a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
+               IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
                   hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
                   hrmean          = MAX(SQRT(Hstar*hi), hi*krdgmin)
@@ -992,5 +908,5 @@
       INTEGER ::   ij                ! horizontal index, combines i and j loops
       INTEGER ::   icells            ! number of cells with aicen > puny
-      REAL(wp) ::   zeps, zindb, zsrdg2   ! local scalar
+      REAL(wp) ::   zindb, zsrdg2   ! local scalar
       REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration
 
@@ -1044,11 +960,13 @@
 
       IF( con_i ) THEN
-         CALL lim_column_sum (jpl,   v_i, vice_init )
-         WRITE(numout,*) ' vice_init  : ', vice_init(jiindx,jjindx)
+         CALL lim_column_sum        (jpl,    v_i,       vice_init )
          CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init )
-         WRITE(numout,*) ' eice_init  : ', eice_init(jiindx,jjindx)
+         DO ji = mi0(jiindx), mi1(jiindx)
+            DO jj = mj0(jjindx), mj1(jjindx)
+               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj)
+               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj)
+            END DO
+         END DO
       ENDIF
-
-      zeps   = 1.e-20_wp
 
       !-------------------------------------------------------------------------------
@@ -1062,5 +980,5 @@
             ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        &
                &                        + opning(ji,jj)                          * rdt_ice
-            IF( ato_i(ji,jj) < -epsi11 ) THEN
+            IF( ato_i(ji,jj) < -epsi10 ) THEN
                neg_ato_i = .TRUE.
             ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error
@@ -1074,9 +992,9 @@
          DO jj = 1, jpj 
             DO ji = 1, jpi
-               IF( ato_i(ji,jj) < -epsi11 ) THEN 
+               IF( ato_i(ji,jj) < -epsi10 ) THEN 
                   WRITE(numout,*) ''  
                   WRITE(numout,*) 'Ridging error: ato_i < 0'
                   WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
-               ENDIF               ! ato_i < -epsi11
+               ENDIF               ! ato_i < -epsi10
             END DO
          END DO
@@ -1120,6 +1038,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( aicen_init(ji,jj,jl1)  >  epsi11  .AND.  athorn(ji,jj,jl1)  >  0._wp       &
-                  .AND. closing_gross(ji,jj) > 0._wp ) THEN
+               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  &
+                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN
                   icells = icells + 1
                   indxi(icells) = ji
@@ -1158,10 +1076,10 @@
             afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
 
-            IF (afrac(ji,jj) > kamax + epsi11) THEN  !riging
+            IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging
                large_afrac = .true.
             ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error
                afrac(ji,jj) = kamax
             ENDIF
-            IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting
+            IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting
                large_afrft = .true.
             ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error
@@ -1222,5 +1140,4 @@
             dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
             dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
-            !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice
             opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
 
@@ -1276,6 +1193,6 @@
 
                ! corrected sea water salinity
-               zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - zeps ) )
-               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), zeps )
+               zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) )
+               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 )
 
                ztmelts          = - tmut * zdummy + rtt
@@ -1312,5 +1229,5 @@
                ji = indxi(ij)
                jj = indxj(ij)
-               IF( afrac(ji,jj) > kamax + epsi11 ) THEN 
+               IF( afrac(ji,jj) > kamax + epsi10 ) THEN 
                   WRITE(numout,*) ''
                   WRITE(numout,*) ' ardg > a_i'
@@ -1324,5 +1241,5 @@
                ji = indxi(ij)
                jj = indxj(ij)
-               IF( afrft(ji,jj) > kamax + epsi11 ) THEN 
+               IF( afrft(ji,jj) > kamax + epsi10 ) THEN 
                   WRITE(numout,*) ''
                   WRITE(numout,*) ' arft > a_i'
@@ -1424,12 +1341,17 @@
          fieldid = ' v_i : limitd_me '
          CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) 
-         WRITE(numout,*) ' vice_init  : ', vice_init(jiindx,jjindx)
-         WRITE(numout,*) ' vice_final : ', vice_final(jiindx,jjindx)
 
          CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final )
          fieldid = ' e_i : limitd_me '
          CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 
-         WRITE(numout,*) ' eice_init  : ', eice_init(jiindx,jjindx)
-         WRITE(numout,*) ' eice_final : ', eice_final(jiindx,jjindx)
+
+         DO ji = mi0(jiindx), mi1(jiindx)
+            DO jj = mj0(jjindx), mj1(jjindx)
+               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj)
+               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
+               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj)
+               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
+            END DO
+         END DO
       ENDIF
       !
@@ -1533,5 +1455,5 @@
 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace
-      
+      REAL(wp)                          ::   zmask_glo
 !!gm      REAL(wp) ::   xtmp      ! temporary variable
       !!-------------------------------------------------------------------
@@ -1545,27 +1467,16 @@
          ! Abort model in case of negative area.
          !-----------------------------------------------------------------
-         IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 .AND. ln_nicep ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF ( a_i(ji,jj,jl) .LT. -epsi11 ) THEN
-                     WRITE (numout,*) ' ALERTE 98 ' 
-                     WRITE (numout,*) ' Negative ice area: ji, jj, jl: ', ji, jj,jl
-                     WRITE (numout,*) ' a_i    ', a_i(ji,jj,jl)
-                  ENDIF
-               END DO
-            END DO
-         ENDIF
-
          icells = 0
-         zmask  = 0._wp
+         zmask(:,:)  = 0._wp
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF(  ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0._wp   ) .OR.   &
-                  & ( a_i(ji,jj,jl) .GT. 0._wp   .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) .OR.   &
-                  & ( v_i(ji,jj,jl)  ==  0._wp   .AND. a_i(ji,jj,jl) .GT. 0._wp   ) .OR.   &
-                  & ( v_i(ji,jj,jl) .GT. 0._wp   .AND. v_i(ji,jj,jl) .LT. 1.e-12  )      )   zmask(ji,jj) = 1._wp
-            END DO
-         END DO
-         IF( ln_nicep )   WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean '
+               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   &
+                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   &
+                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   &
+                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp
+            END DO
+         END DO
+         zmask_glo = glob_sum(zmask)
+         !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean '
 
          !-----------------------------------------------------------------
@@ -1579,5 +1490,5 @@
 !!gm                  xtmp = xtmp * unit_fac
                   ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) )
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )
                END DO
             END DO
@@ -1601,5 +1512,5 @@
 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ???????
 
-               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) )
+               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )
 
                !-----------------------------------------------------------------
@@ -1615,10 +1526,10 @@
 
                ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj)
-               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1 - zmask(ji,jj) )
-               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1 - zmask(ji,jj) )
-               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1 - zmask(ji,jj) )
-               t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1 - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)
-               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1 - zmask(ji,jj) )
-               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )
+               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
+               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
+               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
+               t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)
+               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
+               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
                !
             END DO
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4332)
@@ -45,7 +45,6 @@
    PUBLIC   lim_itd_shiftice
 
-   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values
-   REAL(wp) ::   epsi13 = 1.e-13_wp   !
    REAL(wp) ::   epsi10 = 1.e-10_wp   !
+   REAL(wp) ::   epsi06 = 1.e-6_wp   !
 
    !!----------------------------------------------------------------------
@@ -110,22 +109,5 @@
       CALL lim_var_glo2eqv    ! only for info
 
-      !----------------------------------------------------------------------------------------
-      !  4) Computation of trend terms and get back to old values      
-      !----------------------------------------------------------------------------------------
-
-      !- Trend terms
-      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:) 
-      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:)
-      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)  
-      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 
-      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)
-      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:)
-      d_smv_i_thd(:,:,:) = 0._wp
-      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
-
-      ! diag only (clem)
-      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
-
-      IF(ln_ctl) THEN   ! Control print
+     IF(ln_ctl) THEN   ! Control print
          CALL prt_ctl_info(' ')
          CALL prt_ctl_info(' - Cell values : ')
@@ -183,14 +165,5 @@
       ! -------------------------------
       !
-      !- Recover Old values
-      a_i(:,:,:)   = old_a_i (:,:,:)
-      v_s(:,:,:)   = old_v_s (:,:,:)
-      v_i(:,:,:)   = old_v_i (:,:,:)
-      e_s(:,:,:,:) = old_e_s (:,:,:,:)
-      e_i(:,:,:,:) = old_e_i (:,:,:,:)
-      !?? oa_i(:,:,:)  = old_oa_i(:,:,:)
-      IF( num_sal == 2 )   smv_i(:,:,:) = old_smv_i(:,:,:)
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('limitd_th')
+     IF( nn_timing == 1 )  CALL timing_stop('limitd_th')
    END SUBROUTINE lim_itd_th
    !
@@ -281,9 +254,9 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10))     !0 if no ice and 1 if yes
-               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb
-               zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes
-               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb
-               IF( a_i(ji,jj,jl) > 1e-6 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 
+               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes
+               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb
+               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes
+               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb
+               IF( a_i(ji,jj,jl) > epsi06 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 
             END DO
          END DO
@@ -321,5 +294,5 @@
       ! Tricky trick see limitd_me.F90
       ! will be soon removed, CT
-      ! hi_max(kubnd) = 999.99
+      ! hi_max(kubnd) = 99.
       zhbnew(:,:,:) = 0._wp
 
@@ -329,6 +302,6 @@
             ij = nind_j(ji)
             !
-            IF ( ( zht_i_o(ii,ij,jl)  .GT.epsi10 ) .AND. & 
-               ( zht_i_o(ii,ij,jl+1).GT.epsi10 ) ) THEN
+            IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. & 
+               ( zht_i_o(ii,ij,jl+1) .GT. epsi10 ) ) THEN
                !interpolate between adjacent category growth rates
                zslope = ( zdhice(ii,ij,jl+1)     - zdhice(ii,ij,jl) ) / &
@@ -402,5 +375,7 @@
                zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1)
             ELSE
-               zhbnew(ji,jj,kubnd) = hi_max(kubnd)
+               zhbnew(ji,jj,kubnd) = hi_max(kubnd)  
+               !!? clem bug: since hi_max(jpl)=99, this limit is very high 
+               !!? but I think it is erased in fitline subroutine 
             ENDIF
 
@@ -447,5 +422,5 @@
                      * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
                   a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0
-                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem@useless ?
+                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ?
                ENDIF     ! zetamax > 0
                ! ji, a_i > epsi10
@@ -539,5 +514,5 @@
             a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 
             ht_i(ii,ij,1) = hiclim
-            v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem@useless
+            v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless
          ENDIF
       END DO !ji
@@ -602,13 +577,11 @@
       REAL(wp) ::   zdhr         ! 1 / (hR - hL)
       REAL(wp) ::   zwk1, zwk2   ! temporary variables
-      REAL(wp) ::   zacrith      ! critical minimum concentration in an ice category
-      !!------------------------------------------------------------------
-      !
-      zacrith       = 1.0e-6
+      !!------------------------------------------------------------------
+      !
       !
       DO jj = 1, jpj
          DO ji = 1, jpi
             !
-            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > zacrith   &
+            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   &
                &                        .AND. hice(ji,jj)        > 0._wp     ) THEN
 
@@ -1009,7 +982,10 @@
                   !zdaice(ji,jj,jl)  = a_i(ji,jj,jl)
                   !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)
-                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl)/2
-                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2
+                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp
+                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp
                   ! end TECLIM change 
+                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary
+                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi10 ) / ht_i(ji,jj,jl)  
+                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 )
                ENDIF
             END DO                 ! ji
@@ -1038,26 +1014,40 @@
          zshiftflag = 0
 
+!clem-change
+!         DO jj = 1, jpj
+!            DO ji = 1, jpi
+!               IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   &
+!                  ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
+!                  !
+!                  zshiftflag = 1
+!                  zdonor(ji,jj,jl) = jl + 1
+!                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 
+!                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
+!               ENDIF
+!            END DO                 ! ji
+!         END DO                 ! jj
+!
+!         IF(lk_mpp)   CALL mpp_max( zshiftflag )
+!         
+!         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
+!            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
+!            ! Reset shift parameters
+!            zdonor(:,:,jl) = 0
+!            zdaice(:,:,jl) = 0._wp
+!            zdvice(:,:,jl) = 0._wp
+!         ENDIF
+!clem-change
+
+         ! clem-change begin: why not doing that?
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   &
                   ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
-                  !
-                  zshiftflag = 1
-                  zdonor(ji,jj,jl) = jl + 1
-                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 
-                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
+                  ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10
+                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 
                ENDIF
             END DO                 ! ji
          END DO                 ! jj
-
-         IF(lk_mpp)   CALL mpp_max( zshiftflag )
-         
-         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
-            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
-            ! Reset shift parameters
-            zdonor(:,:,jl) = 0
-            zdaice(:,:,jl) = 0._wp
-            zdvice(:,:,jl) = 0._wp
-         ENDIF
+         ! clem-change end
 
       END DO                    ! jl
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4332)
@@ -50,4 +50,5 @@
    PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
 
+   REAL(wp) ::   epsi10 = 1.e-10_wp   !
    REAL(wp) ::   rzero   = 0._wp   ! constant values
    REAL(wp) ::   rone    = 1._wp   ! constant values
@@ -486,12 +487,4 @@
          CALL lbc_lnk( zs12(:,:), 'F', 1. )
 
-!#if defined key_bdy
-!         ! clem: change zs1, zs2, zs12 at the boundary for each iteration
-!         CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 )
-!         CALL lbc_lnk( zs1 (:,:), 'T', 1. )
-!         CALL lbc_lnk( zs2 (:,:), 'T', 1. )
-!         CALL lbc_lnk( zs12(:,:), 'F', 1. )
-!#endif         
-
          ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
          DO jj = k_j1+1, k_jpj-1
@@ -544,4 +537,8 @@
             CALL agrif_rhg_lim2( jter, nevp, 'U' )
 #endif
+#if defined key_bdy
+         ! clem: change u_ice and v_ice at the boundary for each iteration
+         CALL bdy_ice_lim_dyn( 'U' )
+#endif         
 
 !CDIR NOVERRCHK
@@ -572,4 +569,8 @@
             CALL agrif_rhg_lim2( jter, nevp, 'V' )
 #endif
+#if defined key_bdy
+         ! clem: change u_ice and v_ice at the boundary for each iteration
+         CALL bdy_ice_lim_dyn( 'V' )
+#endif         
 
          ELSE 
@@ -599,6 +600,10 @@
             CALL lbc_lnk( v_ice(:,:), 'V', -1. )
 #if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nevp , 'V' )
-#endif
+            CALL agrif_rhg_lim2( jter, nevp, 'V' )
+#endif
+#if defined key_bdy
+         ! clem: change u_ice and v_ice at the boundary for each iteration
+         CALL bdy_ice_lim_dyn( 'V' )
+#endif         
 
 !CDIR NOVERRCHK
@@ -632,12 +637,11 @@
             CALL agrif_rhg_lim2( jter, nevp, 'U' )
 #endif
+#if defined key_bdy
+         ! clem: change u_ice and v_ice at the boundary for each iteration
+         CALL bdy_ice_lim_dyn( 'U' )
+#endif         
 
          ENDIF
          
-!#if defined key_bdy
-!         ! clem: change u_ice and v_ice at the boundary for each iteration
-!         CALL bdy_ice_lim_dyn( 1 )
-!#endif         
-
          IF(ln_ctl) THEN
             !---  Convergence test.
@@ -666,6 +670,6 @@
 !CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1
-            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
-            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
+            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 
+            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )
             zdummy = vt_i(ji,jj)
             IF ( zdummy .LE. hminrhg ) THEN
@@ -684,11 +688,12 @@
 #if defined key_bdy
       ! clem: change u_ice and v_ice at the boundary
-      CALL bdy_ice_lim_dyn( 1 )
+      CALL bdy_ice_lim_dyn( 'U' )
+      CALL bdy_ice_lim_dyn( 'V' )
 #endif         
 
       DO jj = k_j1+1, k_jpj-1 
          DO ji = fs_2, fs_jpim1
-            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
-            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
+            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 
+            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )
             zdummy = vt_i(ji,jj)
             IF ( zdummy .LE. hminrhg ) THEN
@@ -714,6 +719,6 @@
             !- zdd(:,:), zdt(:,:): divergence and tension at centre 
             !- zds(:,:): shear on northeast corner of grid cells
-            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
-            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
+            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 
+            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )
             zdummy = vt_i(ji,jj)
             IF ( zdummy .LE. hminrhg ) THEN
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4332)
@@ -49,5 +49,4 @@
    PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim
 
-   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
    REAL(wp)  ::   rzero  = 0._wp    
    REAL(wp)  ::   rone   = 1._wp
@@ -145,5 +144,5 @@
 
             !   computation the solar flux at ocean surface
-            IF (lk_cpl) THEN ! be carfeful: not being tested yet
+            IF (lk_cpl) THEN ! be carfeful: not been tested yet
                ! original line
                !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj)
@@ -188,33 +187,4 @@
             qns(ji,jj) = zfcm2 - fdtcn(ji,jj)                        ! non solar heat flux
             !                           ! fdtcn : turbulent oceanic heat flux
-
-            !!gm   this IF prevents the vertorisation of the whole loop
-          !  IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN
-          !     WRITE(numout,*) ' lim_sbc : heat fluxes '
-          !     WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx)
-          !     WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx)
-          !     WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx)
-          !     WRITE(numout,*)
-          !     WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx)
-          !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
-          !     WRITE(numout,*) ' ifral     : ', ifral
-          !     WRITE(numout,*) ' ial       : ', ial  
-          !     WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx)
-          !     WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx)
-          !     !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice
-          !     !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice
-          !     WRITE(numout,*) ' ifrdv     : ', ifrdv
-          !     WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx)
-          !     WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx)
-          !     !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice
-          !     !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice
-          !     WRITE(numout,*) ' '
-          !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
-          !     WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx)
-          !     WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx)
-          !     WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx)
-          !     WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
-          !  ENDIF
-            !!gm   end
          END DO
       END DO
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4332)
@@ -50,9 +50,5 @@
    PUBLIC   lim_thd_init   ! called by iceini module
 
-   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values
-   REAL(wp) ::   epsi16 = 1e-16_wp   !
-   REAL(wp) ::   epsi10 = 1e-10_wp   !
-   REAL(wp) ::   epsi06 = 1e-06_wp   !
-   REAL(wp) ::   epsi04 = 1e-04_wp   !
+   REAL(wp) ::   epsi10 = 1.e-10_wp   !
    REAL(wp) ::   zzero  = 0._wp      !
    REAL(wp) ::   zone   = 1._wp      !
@@ -126,7 +122,7 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i )
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i )
                   !0 if no ice and 1 if yes
-                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) )  ) 
+                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 )  )
                   !convert units ! very important that this line is here
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
@@ -138,7 +134,7 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s )
+                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s )
                   !0 if no ice and 1 if yes
-                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) )  ) 
+                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 )  )
                   !convert units ! very important that this line is here
                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb 
@@ -147,21 +143,4 @@
          END DO
       END DO
-
-      !-----------------------------
-      ! 1.3) Set some dummies to 0
-      !-----------------------------
-      !clem rdvosif(:,:) = 0.e0   ! variation of ice volume at surface
-      !clem rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom
-      !clem fdvolif(:,:) = 0.e0   ! total variation of ice volume
-      !clem rdvonif(:,:) = 0.e0   ! lateral variation of ice volume
-      !clem fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice
-      !clem ffltbif(:,:) = 0.e0   ! linked with fstric
-      !clem qfvbq  (:,:) = 0.e0   ! linked with fstric
-      !clem rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area
-      !clem rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area
-      !clem hicifp (:,:) = 0.e0   ! daily thermodynamic ice production. 
-      !clem sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean
-      !clem fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean
-      !clem sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay
 
       !-----------------------------------
@@ -182,7 +161,5 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-            phicif(ji,jj)  = vt_i(ji,jj)
-            pfrld(ji,jj)   = 1.0 - at_i(ji,jj)
-            zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) )
+            zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) )
             !
             !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
@@ -195,8 +172,8 @@
 
             ! here the drag will depend on ice thickness and type (0.006)
-            fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 
+            fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
             ! also category dependent
             !           !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 
-            qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice
+            qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice
             !                       
             !           !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 
@@ -214,5 +191,5 @@
             zpareff        = 1.0 - zinda * zfntlat
             != 0 if ice and positive heat budget and 1 if one of those two is false
-            zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 )
+            zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) )
             !
             ! Heat budget of the lead, energy transferred from ice to ocean
@@ -221,8 +198,8 @@
             !
             ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx)
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
             !
             ! oceanic heat flux (limthd_dh)
-            fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )
+            fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) )
             !
          END DO
@@ -240,5 +217,5 @@
          ENDIF
 
-         zareamin = 1.e-10
+         zareamin = epsi10
          nbpb = 0
          DO jj = 1, jpj
@@ -248,10 +225,16 @@
                   npb(nbpb) = (jj - 1) * jpi + ji
                ENDIF
-               ! debug point to follow
-               IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN
-                  jiindex_1d = nbpb
-               ENDIF
-            END DO
-         END DO
+            END DO
+         END DO
+
+         ! debug point to follow
+         jiindex_1d = 0
+         IF( ln_nicep ) THEN
+            DO ji = mi0(jiindx), mi1(jiindx)
+               DO jj = mj0(jjindx), mj1(jjindx)
+                  jiindex_1d = (jj - 1) * jpi + ji
+               END DO
+            END DO
+         ENDIF
 
          !------------------------------------------------------------------------------!
@@ -315,6 +298,6 @@
             !--------------------------------
 
-            IF( con_i )   CALL lim_thd_enmelt( 1, nbpb )   ! computes sea ice energy of melting
-            IF( con_i )   CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl )
+            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_enmelt( 1, nbpb )   ! computes sea ice energy of melting
+            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl )
 
             !                                 !---------------------------------!
@@ -324,6 +307,6 @@
             CALL lim_thd_enmelt( 1, nbpb )    ! computes sea ice energy of melting compulsory for limthd_dh
 
-            IF( con_i )   CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 
-            IF( con_i )   CALL lim_thd_con_dif( 1 , nbpb , jl )
+            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 
+            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_con_dif( 1 , nbpb , jl )
 
             !                                 !---------------------------------!
@@ -340,6 +323,6 @@
 
             !           CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
-            IF( con_i )   CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 
-            IF( con_i )   CALL lim_thd_con_dh ( 1 , nbpb , jl )
+            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 
+            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_con_dh ( 1 , nbpb , jl )
 
             !--------------------------------
@@ -431,5 +414,5 @@
 !clem@mv-to-itd    dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
 
-      IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
+      IF( con_i .AND. jiindex_1d > 0 )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
 
       IF(ln_ctl) THEN            ! Control print
@@ -522,8 +505,8 @@
       END DO
       !
-      IF(lwp) WRITE(numout,*) ' lim_thd_glohec '
-      IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice
-      IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice
-      IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice
+      WRITE(numout,*) ' lim_thd_glohec '
+      WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice
+      WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice
+      WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice
       !
    END SUBROUTINE lim_thd_glohec
@@ -611,7 +594,7 @@
       WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit
 
-      IF (jiindex_1D.GT.0) WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d)
-      IF (jiindex_1D.GT.0) WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl)
-      IF (jiindex_1D.GT.0) WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d)
+      WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d)
+      WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl)
+      WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d)
 
       !---------------------------------------
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4332)
@@ -32,9 +32,9 @@
    PUBLIC   lim_thd_dh   ! called by lim_thd
 
-   REAL(wp) ::   epsi20 = 1e-20   ! constant values
-   REAL(wp) ::   epsi13 = 1e-13   !
-   REAL(wp) ::   epsi16 = 1e-16   !
-   REAL(wp) ::   zzero  = 0.e0    !
-   REAL(wp) ::   zone   = 1.e0    !
+   REAL(wp) ::   epsi20 = 1.e-20   ! constant values
+   REAL(wp) ::   epsi10 = 1.e-10   !
+   REAL(wp) ::   epsi13 = 1.e-13   !
+   REAL(wp) ::   zzero  = 0._wp    !
+   REAL(wp) ::   zone   = 1._wp    !
 
    !!----------------------------------------------------------------------
@@ -125,4 +125,5 @@
       INTEGER  ::   num_iter_max, numce_dh
       REAL(wp) ::   meance_dh
+      REAL(wp) ::   zinda 
       REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt
       REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i
@@ -288,7 +289,7 @@
       END DO
 
-      !                     !-------------------
-      IF( con_i ) THEN      ! Conservation test
-         !                  !-------------------
+      !                                          !-------------------
+      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test
+         !                                       !-------------------
          numce_dh  = 0
          meance_dh = 0._wp
@@ -494,7 +495,7 @@
       END DO ! jk
 
-      !                     !-------------------
-      IF( con_i ) THEN      ! Conservation test
-      !                     !-------------------
+      !                                          !-------------------
+      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test
+      !                                          !-------------------
          DO ji = kideb, kiut
             IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN
@@ -544,5 +545,6 @@
          sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice
          !                     ! excessive energy is sent to lateral ablation
-         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice
+         zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) )
+         fsup(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice
       END DO
 
@@ -577,5 +579,6 @@
 
          ! energy of melting of remaining snow
-         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 )
+         zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) )
+         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda
          zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 )
          zhnfi     =  zhni + zdhnm
@@ -613,5 +616,5 @@
          ENDIF
          !
-         ! diagnostic ( bottom ice growth )
+         ! diagnostic 
          ii = MOD( npb(ji) - 1, jpi ) + 1
          ij = ( npb(ji) - 1 ) / jpi + 1
@@ -689,10 +692,10 @@
          ! clem: new salinity difference stored (to be used in limthd_ent.F90)
          IF (  num_sal == 2  ) THEN
-            i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) )
+            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) )
             ! salinity dif due to snow-ice formation
-            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch     
+            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch     
             ! salinity dif due to bottom growth 
             IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN
-               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch
+               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch
             ENDIF
          ENDIF
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4332)
@@ -31,7 +31,5 @@
    PUBLIC   lim_thd_dif   ! called by lim_thd
 
-   REAL(wp) ::   epsi20 = 1e-20     ! constant values
-   REAL(wp) ::   epsi13 = 1e-13     ! constant values
-
+   REAL(wp) ::   epsi10      =  1.e-10_wp    !
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
@@ -107,5 +105,4 @@
       INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation
       INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow
-      REAL(wp) ::   zeps      =  1.e-10_wp    !
       REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
       REAL(wp) ::   zg1       =  2._wp        !
@@ -114,5 +111,4 @@
       REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow
       REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
-      REAL(wp) ::   zht_smin  =  1.e-4_wp     ! minimum snow depth
       REAL(wp) ::   ztmelt_i    ! ice melting temperature
       REAL(wp) ::   zerritmax   ! current maximal error on temperature 
@@ -312,5 +308,5 @@
          IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula
             DO ji = kideb , kiut
-               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-zeps,t_i_b(ji,1)-rtt)
+               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)
                ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin)
             END DO
@@ -318,5 +314,5 @@
                DO ji = kideb , kiut
                   ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  &
-                     MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
+                     MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
                   ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
                END DO
@@ -326,5 +322,5 @@
          IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
             DO ji = kideb , kiut
-               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -zeps, t_i_b(ji,1)-rtt )   &
+               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   &
                   &                   - 0.011_wp * ( t_i_b(ji,1) - rtt )  
                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
@@ -333,5 +329,5 @@
                DO ji = kideb , kiut
                   ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
-                     &                                  / MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   &
+                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   &
                      &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )  
                   ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
@@ -339,5 +335,5 @@
             END DO
             DO ji = kideb , kiut
-               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-zeps,t_bo_b(ji)-rtt)   &
+               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   &
                   &                        - 0.011_wp * ( t_bo_b(ji) - rtt )  
                ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
@@ -352,6 +348,6 @@
 
             !-- Snow kappa factors
-            zkappa_s(ji,0)         = rcdsn / MAX(zeps,zh_s(ji))
-            zkappa_s(ji,nlay_s)    = rcdsn / MAX(zeps,zh_s(ji))
+            zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji))
+            zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji))
          END DO
 
@@ -359,5 +355,5 @@
             DO ji = kideb , kiut
                zkappa_s(ji,layer)  = 2.0 * rcdsn / &
-                  MAX(zeps,2.0*zh_s(ji))
+                  MAX(epsi10,2.0*zh_s(ji))
             END DO
          END DO
@@ -367,13 +363,13 @@
                !-- Ice kappa factors
                zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ &
-                  MAX(zeps,2.0*zh_i(ji)) 
-            END DO
-         END DO
-
-         DO ji = kideb , kiut
-            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(zeps,zh_i(ji))
-            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji))
+                  MAX(epsi10,2.0*zh_i(ji)) 
+            END DO
+         END DO
+
+         DO ji = kideb , kiut
+            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))
+            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))
             !-- Interface
-            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, &
+            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, &
                (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji)))
             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) &
@@ -389,7 +385,7 @@
                ztitemp(ji,layer)   = t_i_b(ji,layer)
                zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
-                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps)
+                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)
                zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &
-                  zeps)
+                  epsi10)
             END DO
          END DO
@@ -398,5 +394,5 @@
             DO ji = kideb , kiut
                ztstemp(ji,layer) = t_s_b(ji,layer)
-               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),zeps)
+               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)
             END DO
          END DO
@@ -707,5 +703,5 @@
       END DO  ! End of the do while iterative procedure
 
-      IF( ln_nicep ) THEN
+      IF( ln_nicep .AND. lwp ) THEN
          WRITE(numout,*) ' zerritmax : ', zerritmax
          WRITE(numout,*) ' nconv     : ', nconv
@@ -732,5 +728,5 @@
       ! Heat conservation       !
       !-------------------------!
-      IF( con_i ) THEN
+      IF( con_i .AND. jiindex_1d > 0 ) THEN
          DO ji = kideb, kiut
             ! Upper snow value
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4332)
@@ -36,8 +36,6 @@
    PUBLIC   lim_thd_ent         ! called by lim_thd
 
-   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values
-   REAL(wp) ::   epsi13 = 1e-13_wp   !
-   REAL(wp) ::   epsi10 = 1e-10_wp   !
-   REAL(wp) ::   epsi06 = 1e-06_wp   !
+   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values
+   REAL(wp) ::   epsi10 = 1.e-10_wp   !
    REAL(wp) ::   zzero  = 0._wp      !
    REAL(wp) ::   zone   = 1._wp      !
@@ -122,4 +120,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zhl0      ! old and new layer thicknesses 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zrl01
+
+      REAL(wp) ::   zinda 
       !!-------------------------------------------------------------------
 
@@ -391,5 +391,6 @@
          DO layer1 = ntop1, nbot1
             DO ji = kideb, kiut
-               zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   &
+               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) )
+               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   &
                   &                 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 
                q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0)   &
@@ -407,5 +408,5 @@
       END DO
 
-      IF ( con_i ) THEN
+      IF ( con_i .AND. jiindex_1d > 0 ) THEN
          DO ji = kideb, kiut
             IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN
@@ -429,5 +430,6 @@
       DO jk = 1, nlay_s
          DO ji = kideb, kiut
-            q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , epsi20 )
+            zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) )        
+            q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 )
          END DO !ji
       END DO !jk  
@@ -555,5 +557,5 @@
          DO ji = kideb, kiut
             ! Heat conservation
-            zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) &
+            zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) &
                &                                   * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) )
          END DO
@@ -601,9 +603,10 @@
          DO layer1 = ntop1, nbot1
             DO ji = kideb, kiut
-               zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) &
+               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) )
+               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) &
                   - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10))
                q_i_b(ji,layer1) = q_i_b(ji,layer1) & 
                   + zrl01(layer1,layer0)*qm0(ji,layer0) &
-                  * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) &
+                  * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) &
                   * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0)))
             END DO
@@ -621,5 +624,5 @@
       END DO
       !
-      IF ( con_i ) THEN
+      IF ( con_i .AND. jiindex_1d > 0 ) THEN
          DO ji = kideb, kiut
             IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN
@@ -646,5 +649,6 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut
-            q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , epsi20 )
+            zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) )
+            q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 )
          END DO !ji
       END DO !jk  
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4332)
@@ -36,10 +36,5 @@
    PUBLIC lim_thd_lac     ! called by lim_thd
 
-   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values
-   REAL(wp) ::   epsi13 = 1e-13_wp   !
-   REAL(wp) ::   epsi11 = 1e-11_wp   !
-   REAL(wp) ::   epsi10 = 1e-10_wp   !
-   REAL(wp) ::   epsi06 = 1e-06_wp   !
-   REAL(wp) ::   epsi03 = 1e-03_wp   !
+   REAL(wp) ::   epsi10 = 1.e-10_wp   !
    REAL(wp) ::   zzero  = 0._wp      !
    REAL(wp) ::   zone   = 1._wp      !
@@ -160,5 +155,5 @@
                   !Energy of melting q(S,T) [J.m-3]
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i )
-                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes
+                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
                END DO
@@ -286,10 +281,21 @@
                nbpac = nbpac + 1
                npac( nbpac ) = (jj - 1) * jpi + ji
-               IF( ji == jiindx .AND. jj == jjindx )   jiindex_1d = nbpac
             ENDIF
          END DO
       END DO
 
-      IF( ln_nicep )   WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
+      ! debug point to follow
+      jiindex_1d = 0
+      IF( ln_nicep ) THEN
+         DO ji = mi0(jiindx), mi1(jiindx)
+            DO jj = mj0(jjindx), mj1(jjindx)
+               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN
+                  jiindex_1d = (jj - 1) * jpi + ji
+               ENDIF
+            END DO
+         END DO
+      ENDIF
+   
+      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
 
       !------------------------------
@@ -499,13 +505,13 @@
          END DO
 
-         IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
+         IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
             DO ji = 1, nbpac
                zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
-               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem
-               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 )
-            END DO
-         END DO
-         IF( ln_nicep )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
+               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) )  ! clem
+               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 )
+            END DO
+         END DO
+         IF( ln_nicep .AND. jiindex_1d > 0 )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
 
          !---------------------------------
@@ -649,8 +655,12 @@
          !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 
          IF( ln_nicep ) THEN
-            WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindx,jjindx)
-            WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindx,jjindx)
-            WRITE(numout,*) ' et_i_init : ', et_i_init(jiindx,jjindx)
-            WRITE(numout,*) ' et_i_final: ', et_i_final(jiindx,jjindx)
+            DO ji = mi0(jiindx), mi1(jiindx)
+               DO jj = mj0(jjindx), mj1(jjindx)
+                  WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
+                  WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
+                  WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
+                  WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
+               END DO
+            END DO
          ENDIF
          !
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4332)
@@ -36,13 +36,7 @@
    PUBLIC   lim_trp    ! called by ice_step
 
-   REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values
-   REAL(wp)  ::   epsi03 = 1.e-03_wp  
    REAL(wp)  ::   epsi10 = 1.e-10_wp  
-   REAL(wp)  ::   epsi16 = 1.e-16_wp
-   REAL(wp)  ::   epsi20 = 1.e-20_wp
    REAL(wp)  ::   rzero  = 0._wp   
    REAL(wp)  ::   rone   = 1._wp
-
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) ::   zs0e
 
    !! * Substitution
@@ -448,6 +442,6 @@
 
                   ! Switches and dummy variables
-                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
-                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
+                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi10 )
+                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi10 )
                   zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
                   zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
@@ -460,5 +454,5 @@
                   ENDIF
 
-                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp  ) * a_i(ji,jj,jl)
+                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp  ) * a_i(ji,jj,jl)
                   oa_i (ji,jj,jl)  = zindic * zage 
 
@@ -500,6 +494,6 @@
                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
                !
-               zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi16 ) )
-               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness
+               zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) )
+               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness
             END DO
          END DO
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4332)
@@ -48,10 +48,5 @@
    PUBLIC   lim_update1   ! routine called by ice_step
 
-      REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants
-      REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       -
-      REAL(wp)  ::   epsi03 = 1.e-03_wp   !    -       -
       REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
-      REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       -
-      REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       -
       REAL(wp)  ::   rzero  = 0._wp       !    -       -
       REAL(wp)  ::   rone   = 1._wp       !    -       -
@@ -108,37 +103,17 @@
       !------------------------------------------------------------------------------
 
-      !---------------------
-      ! Ice dynamics  
-      !---------------------
-      u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:)
-      v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:)
- 
-      !-----------------------------
-      ! Update ice and snow volumes  
-      !-----------------------------
-      DO jl = 1, jpl
-         v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_trp(:,:,jl) 
-         v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_trp(:,:,jl)
-      END DO
-
-
-      !---------------------------------------------
-      ! Ice concentration and ice heat content
-      !---------------------------------------------
-      a_i (:,:,:)  = a_i (:,:,:) + d_a_i_trp(:,:,:)
-      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:)  
-
-      !------------------------------
-      ! Snow temperature and ice age
-      !------------------------------
-      e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:)
-      oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_trp(:,:,:)  
-
-      !--------------
-      ! Ice salinity    
-      !--------------
-      IF(  num_sal == 2  ) THEN      ! general case
-         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_trp(:,:,:)
-      ENDIF
+      !-----------------
+      !  Trend terms
+      !-----------------
+      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:)
+      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:)
+      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:)
+      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)  
+      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)   
+      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)  
+      d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:)
+      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:)
+      d_smv_i_trp(:,:,:)   = 0._wp
+      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:)
 
       ! mass and salt flux init (clem)
@@ -160,94 +135,52 @@
       CALL lim_var_glo2eqv
 
-      !---------------------------------
-      ! Classify the pathological cases
-      !---------------------------------
-      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
-      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
-      ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               patho_case(ji,jj,jl) = 1
-               IF( v_i(ji,jj,jl) .LT. 0.0 ) THEN
-                  patho_case(ji,jj,jl) = 3
-               ENDIF
-               IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
-                  patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free cell
-               ENDIF
-            END DO
-         END DO
-      END DO
-
       !--------------------------------------
       ! 2. Review of all pathological cases
       !--------------------------------------
 
+! clem: useless now
       !-------------------------------------------
       ! 2.1) Advection of ice in an ice-free cell
       !-------------------------------------------
       ! should be removed since it is treated after dynamics now
-
-!      !IF( ln_nicep ) THEN  
-!         WRITE(numout,*) ' limupdate1 - before h correction '
-!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' at_i : ', at_i(12,44)
-!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
-!      !ENDIF
+!      zhimax = 5._wp
+!      ! first category
+!      DO jj = 1, jpj
+!         DO ji = 1, jpi
+!            !--- the thickness of such an ice is often out of bounds
+!            !--- thus we recompute a new area while conserving ice volume
+!            zat_i_old = SUM( old_a_i(ji,jj,:) )
+!            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) ) 
+!            IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) &
+!              &   .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) &
+!              &   .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
+!               ht_i(ji,jj,1) = hi_max(1) * 0.5_wp
+!               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
+!            ENDIF
+!         END DO
+!      END DO
 !
-      zhimax = .3_wp
-      ! first category
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !--- the thickness of such an ice is often out of bounds
-            !--- thus we recompute a new area while conserving ice volume
-            zat_i_old = SUM( old_a_i(ji,jj,:) )
-            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) ) 
-            IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) &
-                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) &
-                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
-               ht_i(ji,jj,1) = hi_max(1) / 2.0
-               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
-            ENDIF
-         END DO
-      END DO
-
-
-!      !IF( ln_nicep ) THEN  
-!         at_i(:,:) = 0._wp
-!         DO jl = 1, jpl
-!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-!         END DO
-!         WRITE(numout,*) ' limupdate1 - after h correction 1 '
-!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' at_i : ', at_i(12,44)
-!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
-!      !ENDIF
-
-      zhimax = 1._wp
-      ! other categories
-      DO jl = 2, jpl
-         jm = ice_types(jl)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 
-               ! this correction is very tricky... sometimes, advection gets wrong i don't know why
-               ! it makes problems when the advected volume and concentration do not seem to be 
-               ! related with each other
-               ! the new thickness is sometimes very big!
-               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
-               ! which of course is plausible
-               ! but fuck! it fucks everything up :)
-               IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) &
-                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN
-                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0
-                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
-               ENDIF
-            END DO ! ji
-         END DO !jj
-      END DO !jl
+!      zhimax = 20._wp
+!      ! other categories
+!      DO jl = 2, jpl
+!         jm = ice_types(jl)
+!         DO jj = 1, jpj
+!            DO ji = 1, jpi
+!               zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) ) 
+!               ! this correction is very tricky... sometimes, advection gets wrong i don't know why
+!               ! it makes problems when the advected volume and concentration do not seem to be 
+!               ! related with each other
+!               ! the new thickness is sometimes very big!
+!               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
+!               ! which of course is plausible
+!               ! but fuck! it fucks everything up :)
+!               IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) &
+!                  &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN
+!                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp
+!                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
+!               ENDIF
+!            END DO ! ji
+!         END DO !jj
+!      END DO !jl
      
       at_i(:,:) = 0._wp
@@ -255,12 +188,4 @@
          at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
       END DO
-
-!      !IF( ln_nicep ) THEN  
-!         WRITE(numout,*) ' limupdate1 - after h correction 2 '
-!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' at_i : ', at_i(12,44)
-!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
-!      !ENDIF
 
       !----------------------------------------------------
@@ -285,8 +210,8 @@
 
                !switches
-               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
+               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 
                !switch = 1 if a_i > 1e-06 and 0 if not
-               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
-               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
+               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not
+               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not
                ! bug fix 25 avril 2007
                zindb         = zindb*zindic
@@ -332,5 +257,4 @@
                !--------------------------------------------
                ! if greater than 0, ice concentration cannot be smaller than 1e-10
-               !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
                a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl)
 
@@ -352,5 +276,5 @@
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
+                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
                   e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) )
                END DO ! ji
@@ -370,11 +294,11 @@
          DO ji = 1, jpi
             z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )        
-            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 
+            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 
             DO jl  = 1, jpl
-               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
+               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb
                a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i )
                !
-               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
-               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda
+               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 
+               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda
                !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used
             END DO
@@ -412,6 +336,6 @@
                      !clem smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) )
                      smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
-                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) )
-                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
+                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) )
+                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
                   END DO ! ji
                END DO ! jj
@@ -512,5 +436,4 @@
             CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ')
             CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ')
-            CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update1  : Path. case  : ')
 
             DO jk = 1, nlay_i
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4332)
@@ -45,9 +45,5 @@
    PUBLIC   lim_update2   ! routine called by ice_step
 
-      REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants
-      REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       -
       REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
-      REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       -
-      REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       -
       REAL(wp)  ::   rzero  = 0._wp       !    -       -
       REAL(wp)  ::   rone   = 1._wp       !    -       -
@@ -102,34 +98,18 @@
       CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem
 
-      !------------------------------------------------------------------------------
-      ! 1. Update of Global variables                                               |
-      !------------------------------------------------------------------------------
-
-      !-----------------------------
-      ! Update ice and snow volumes  
-      !-----------------------------
-      DO jl = 1, jpl
-         v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_thd(:,:,jl) 
-         v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_thd(:,:,jl)
-      END DO
- 
-      !---------------------------------------------
-      ! Ice concentration and ice heat content
-      !---------------------------------------------
-      a_i (:,:,:) = a_i (:,:,:) + d_a_i_thd(:,:,:)
-      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:)  
- 
-      !------------------------------
-      ! Snow temperature and ice age
-      !------------------------------
-      e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_thd (:,:,:,:)
-      oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_thd(:,:,:)
-
-      !--------------
-      ! Ice salinity    
-      !--------------
-      IF(  num_sal == 2  ) THEN      ! general case
-         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:)
-      ENDIF
+      !----------------------------------------------------------------------------------------
+      ! 1. Computation of trend terms      
+      !----------------------------------------------------------------------------------------
+      !- Trend terms
+      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:) 
+      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:)
+      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)  
+      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 
+      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)
+      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:)
+      d_smv_i_thd(:,:,:) = 0._wp
+      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
+      ! diag only (clem)
+      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
 
       ! mass and salt flux init (clem)
@@ -151,98 +131,52 @@
       CALL lim_var_glo2eqv
 
-      !---------------------------------
-      ! Classify the pathological cases
-      !---------------------------------
-      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
-      ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)
-      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
-      ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation 
-      ! with negative advection, very pathological )
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               patho_case(ji,jj,jl) = 1
-               IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN
-                  IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN 
-                     patho_case(ji,jj,jl) = 2
-                  ENDIF
-               ELSE
-                  patho_case(ji,jj,jl) = 3
-                  IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN 
-                     patho_case(ji,jj,jl) = 4
-                  ENDIF
-               ENDIF
-             END DO
-         END DO
-      END DO
-
       !--------------------------------------
       ! 2. Review of all pathological cases
       !--------------------------------------
+
+! clem: useless now
       !-------------------------------------------
       ! 2.1) Advection of ice in an ice-free cell
       !-------------------------------------------
       ! should be removed since it is treated after dynamics now
-
-!      !IF( ln_nicep ) THEN  
-!         WRITE(numout,*) ' limupdate2 - before h correction '
-!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' at_i : ', at_i(12,44)
-!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
-!      !ENDIF
-
-      zhimax = .3_wp
-      ! first category
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !--- the thickness of such an ice is often out of bounds
-            !--- thus we recompute a new area while conserving ice volume
-            zat_i_old = SUM( old_a_i(ji,jj,:) )
-            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1)) - epsi10 ) ) 
-            IF ( ( ABS(d_v_i_thd(ji,jj,1))/MAX(ABS(d_a_i_thd(ji,jj,1)),epsi10)*zindb .GT. zhimax) &
-                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) &
-                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
-               ht_i(ji,jj,1) = hi_max(1) / 2.0
-               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
-            ENDIF
-         END DO
-      END DO
-
-!      !IF( ln_nicep ) THEN  
-!         at_i(:,:) = 0._wp
-!         DO jl = 1, jpl
-!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+!      zhimax = 5._wp
+!      ! first category
+!      DO jj = 1, jpj
+!         DO ji = 1, jpi
+!            !--- the thickness of such an ice is often out of bounds
+!            !--- thus we recompute a new area while conserving ice volume
+!            zat_i_old = SUM( old_a_i(ji,jj,:) )
+!            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1) ) - epsi10 ) ) 
+!            IF ( ( ABS( d_v_i_thd(ji,jj,1) ) / MAX( ABS( d_a_i_thd(ji,jj,1) ),epsi10 ) * zindb .GT. zhimax ) &
+!               &  .AND. ( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) &
+!               &  .AND. ( zat_i_old .LT. 1.e-6 ) )  THEN ! new line
+!               ht_i(ji,jj,1) = hi_max(1) * 0.5_wp
+!               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)
+!            ENDIF
 !         END DO
-!         WRITE(numout,*) ' limupdate2 - after h correction 1 '
-!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' at_i : ', at_i(12,44)
-!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
-!      !ENDIF
-! 
-      zhimax = 1._wp
-      ! other categories
-      DO jl = 2, jpl
-         jm = ice_types(jl)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 
-              ! this correction is very tricky... sometimes, advection gets wrong i don't know why
-               ! it makes problems when the advected volume and concentration do not seem to be 
-               ! related with each other
-               ! the new thickness is sometimes very big!
-               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
-               ! which of course is plausible
-               ! but fuck! it fucks everything up :)
-               IF ( (ABS(d_v_i_thd(ji,jj,jl))/MAX(ABS(d_a_i_thd(ji,jj,jl)),epsi10)*zindb .GT. zhimax) &
-                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN
-                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0
-                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
-               ENDIF
-            END DO ! ji
-         END DO !jj
-      END DO !jl
+!      END DO
+
+!      zhimax = 20._wp
+!      ! other categories
+!      DO jl = 2, jpl
+!         jm = ice_types(jl)
+!         DO jj = 1, jpj
+!            DO ji = 1, jpi
+!               zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 
+!              ! this correction is very tricky... sometimes, advection gets wrong i don't know why
+!               ! it makes problems when the advected volume and concentration do not seem to be 
+!               ! related with each other
+!               ! the new thickness is sometimes very big!
+!               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
+!               ! which of course is plausible
+!               ! but fuck! it fucks everything up :)
+!               IF ( ( ABS( d_v_i_thd(ji,jj,jl) ) / MAX( ABS( d_a_i_thd(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) &
+!                  &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN
+!                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp
+!                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)
+!               ENDIF
+!            END DO ! ji
+!         END DO !jj
+!      END DO !jl
      
       at_i(:,:) = 0._wp
@@ -250,12 +184,4 @@
          at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
       END DO
-
-!      !IF( ln_nicep ) THEN  
-!         WRITE(numout,*) ' limupdate2 - after h correction 2 '
-!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' at_i : ', at_i(12,44)
-!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl)
-!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)
-!      !ENDIF
 
       !----------------------------------------------------
@@ -304,7 +230,4 @@
                      zesum = zesum + e_i(ji,jj,jk,jl)
                   END DO
-                  !clem IF (ind_im .LT. nlay_i ) THEN
-                  !clem   smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) )
-                  !clem ENDIF
                   ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i)
                   v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
@@ -368,5 +291,6 @@
             DO ji = 1, jpi
                ! snow energy of melting
-               ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), 1.0e-6 )  ! snow energy of melting
+               zinda   =  MAX( 0._wp, SIGN( 1._wp, v_s(ji,jj,jl) - epsi10 ) )
+               ze_s = zinda * e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), epsi10 )  ! snow energy of melting
 
                ! If snow energy of melting smaller then Lf
@@ -401,8 +325,8 @@
 
                !switches
-               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
+               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 
                !switch = 1 if a_i > 1e-06 and 0 if not
-               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
-               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
+               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not
+               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not
                ! bug fix 25 avril 2007
                zindb         = zindb*zindic
@@ -444,5 +368,4 @@
                !--- 2.7 Correction to ice concentrations 
                !--------------------------------------------
-               !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
                a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl)
 
@@ -464,5 +387,5 @@
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) 
+                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
                   e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) )
                END DO ! ji
@@ -478,7 +401,7 @@
                !--- 2.12 Constrain the thickness of the smallest category above 5 cm
                !----------------------------------------------------------------------
-               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
-               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06)
-               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi20 ) )
+               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 
+               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi10)
+               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi10 ) )
                ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh
                ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh
@@ -502,11 +425,11 @@
          DO ji = 1, jpi
             z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )        
-            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 
+            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 
             DO jl  = 1, jpl
-               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb
+               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb
                a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i )
                !
-               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
-               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda
+               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 
+               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda
                !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used
             END DO
@@ -542,6 +465,6 @@
                      !clem smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) )
                      smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
-                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) )
-                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
+                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) )
+                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
                   END DO ! ji
                END DO ! jj
@@ -663,5 +586,4 @@
             CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update2  : old_oa_i    : ')
             CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ')
-            CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update2  : Path. case  : ')
 
             DO jk = 1, nlay_i
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4332)
@@ -66,9 +66,5 @@
    PUBLIC   lim_var_salprof1d    !
 
-   REAL(wp) ::   epsi20 = 1.e-20_wp   ! module constants
-   REAL(wp) ::   epsi16 = 1.e-16_wp   !    -       -
-   REAL(wp) ::   epsi13 = 1.e-13_wp   !    -       -
    REAL(wp) ::   epsi10 = 1.e-10_wp   !    -       -
-   REAL(wp) ::   epsi06 = 1.e-06_wp   !    -       -
    REAL(wp) ::   zzero = 0.e0        !    -       -
    REAL(wp) ::   zone  = 1.e0        !    -       -
@@ -117,6 +113,6 @@
                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
                !
-               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) ) 
-               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness
+               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) ) 
+               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness
             END DO
          END DO
@@ -138,9 +134,9 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi16 ) ) 
-                  zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) ) 
+                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi10 ) ) 
+                  zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) ) 
                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content
-                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi16 ) * zinda   ! ice salinity
-                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi16 ) * zindb   ! ice age
+                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda   ! ice salinity
+                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * zindb   ! ice age
                END DO
             END DO
@@ -209,6 +205,6 @@
                DO ji = 1, jpi
                   !                                                              ! Energy of melting q(S,T) [J.m-3]
-                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi06 ) * REAL(nlay_i,wp) 
-                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes
+                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
+                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes
                   zq_i    = zq_i * unit_fac * zindb                              !convert units
                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature
@@ -235,6 +231,6 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL(nlay_s,wp)
-                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes
+                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)
+                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes
                   zq_s  = zq_s * unit_fac * zindb                                    ! convert units
                   !
@@ -341,4 +337,5 @@
                   zind01 = ( 1._wp - zind0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) ) 
                   ! If 2.sm_i GE sss_m then zindbal = 1
+                  ! this is to force a constant salinity profile in the Baltic Sea
                   zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
                   zalpha(ji,jj,jl) = zind0  + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )
@@ -434,9 +431,9 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi16 ) )  )
-                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi16 ) )  )
-                  zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi16 )   &
+                  zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  )
+                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
+                  zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   &
                      &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp)
-                  bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi16 )
+                  bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi10 )
                END DO
             END DO
@@ -498,4 +495,5 @@
                zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) ) 
                ! if 2.sm_i GE sss_m then zindbal = 1
+               ! this is to force a constant salinity profile in the Baltic Sea
                zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) )
                !
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4220)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4332)
@@ -52,5 +52,5 @@
    INTEGER            , DIMENSION(jpnoumax) ::   nc  , nca     ! switch for saving field ( = 1 ) or not ( = 0 )
 
-   REAL(wp)  ::   epsi06 = 1e-6_wp
+   REAL(wp)  ::   epsi06 = 1.e-6_wp
    REAL(wp)  ::   zzero  = 0._wp
    REAL(wp)  ::   zone   = 1._wp      
@@ -192,10 +192,10 @@
 
       !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-      IF( ln_nicep ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit
-         WRITE(numout,*) '~~~~~~~ '
-         WRITE(numout,*) ' kindic = ', kindic
-      ENDIF
+      !IF( ln_nicep ) THEN
+      !   WRITE(numout,*)
+      !   WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit
+      !   WRITE(numout,*) '~~~~~~~ '
+      !   WRITE(numout,*) ' kindic = ', kindic
+      !ENDIF
       !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
@@ -564,4 +564,15 @@
       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"      , 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, "iisfxthd", "Salt flux from thermo"   , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iisfxmec", "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, "iice_itd", "Ice concentration by cat", "%"      , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) 
@@ -586,4 +597,15 @@
       CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
 
+      CALL histwrite( kid, "iicebopr", kt, diag_bot_gr        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicedypr", kt, diag_dyn_gr        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicelapr", kt, diag_lat_gr        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesipr", kt, diag_sni_gr        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicerepr", kt, diag_res_pr        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicebome", kt, diag_bot_me        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iicesume", kt, diag_sur_me        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iisfxthd", kt, sfx_thd        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iisfxmec", kt, sfx_mec        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) )
+
       !CALL histwrite( kid, "iice_itd", kt, a_i  , jpi*jpj*jpl, (/1/)  )   ! area
       !CALL histwrite( kid, "iice_hid", kt, ht_i , jpi*jpj*jpl, (/1/)  )   ! thickness
