Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/cfg.txt	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/cfg.txt	(revision 4333)
@@ -1,5 +1,4 @@
 GYRE_BFM OPA_SRC TOP_SRC
 GYRE_PISCES OPA_SRC TOP_SRC
-ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC
 ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
 GYRE OPA_SRC
@@ -10,2 +9,3 @@
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
+ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4333)
@@ -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
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4333)
@@ -90,5 +90,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
@@ -105,8 +105,8 @@
       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(:,:,:)
+      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      
@@ -143,8 +143,8 @@
       WRITE ( numoni, 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
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90	(revision 4333)
@@ -45,5 +45,5 @@
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
       REAL(wp) :: epsi06 = 1.0e-6
-      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
+      REAL(wp) :: zc1, zc2, zc3, zx1, zdh   ! local scalars
       REAL(wp), DIMENSION(0:jpl)   ::   zhi_max         !:Boundary of ice thickness categories in thickness space
  
@@ -85,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
@@ -108,5 +108,4 @@
                zht_i(ji,1) = zhti(ji)
                za_i (ji,1) = zai (ji)
-               
                ! *** case ice is thicker: fill categories >1
             ELSE
@@ -178,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: anches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90	(revision 4332)
+++ 	(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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 4333)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4333)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4333)
@@ -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                             ! 
@@ -261,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
@@ -281,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
@@ -312,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
@@ -354,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
@@ -373,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)
@@ -396,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
@@ -415,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
@@ -487,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.
       !
@@ -601,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)
                      !----------------------------
@@ -620,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
@@ -677,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)              &
@@ -717,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
@@ -790,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
@@ -800,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
@@ -883,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
@@ -934,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)
@@ -988,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
 
@@ -1040,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
 
       !-------------------------------------------------------------------------------
@@ -1058,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
@@ -1070,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
@@ -1116,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
@@ -1154,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
@@ -1218,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
 
@@ -1272,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
@@ -1308,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'
@@ -1320,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'
@@ -1420,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
       !
@@ -1536,5 +1462,5 @@
 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace
-      
+      REAL(wp)                          ::   zmask_glo
 !!gm      REAL(wp) ::   xtmp      ! temporary variable
       !!-------------------------------------------------------------------
@@ -1548,27 +1474,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 '
 
          !-----------------------------------------------------------------
@@ -1582,5 +1497,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
@@ -1604,5 +1519,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) )
 
                !-----------------------------------------------------------------
@@ -1618,10 +1533,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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4333)
@@ -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
@@ -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) ) / &
@@ -449,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
@@ -541,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
@@ -604,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
 
@@ -1015,5 +986,5 @@
                   ! 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) ) / ht_i(ji,jj,jl)  
+                  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
@@ -1043,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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4333)
@@ -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
@@ -426,17 +427,17 @@
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
-               !   &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
-               !   &          / ( 1._wp + alphaevp * dtotel )
-
-               !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
-               !              zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
-               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel )
+               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
+                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
+                  &          / ( 1._wp + alphaevp * dtotel )
+
+               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
+                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )
 
                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
-               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
-                  &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
-               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
-                  &         / ( 1._wp + dtotel )
+               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
+               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
+               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
+               !   &         / ( 1._wp + dtotel )
 
             END DO
@@ -472,12 +473,12 @@
 
                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
-               !   &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
-               !   &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
+               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
+                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
 
                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
-               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
-                  &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
-                  &          / ( 1.0 + dtotel ) 
+               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
+               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+               !   &          / ( 1.0 + dtotel ) 
 
             END DO ! ji
@@ -485,12 +486,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)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4333)
@@ -51,5 +51,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
@@ -147,5 +146,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)
@@ -190,33 +189,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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4333)
@@ -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) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj) * ( 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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4333)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4333)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4333)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4333)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4333)
@@ -36,9 +36,5 @@
    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
@@ -446,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 ) )
@@ -458,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 
 
@@ -498,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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4333)
@@ -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 = 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
-
-
-!      !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 = 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
+!      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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4333)
@@ -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 = 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
-      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 = 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
+!      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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4333)
@@ -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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4333)
@@ -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
       !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
@@ -570,4 +570,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 ) 
@@ -592,4 +603,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
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4333)
@@ -104,7 +104,7 @@
    REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points
 
-#if defined key_lim2
-   CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim2      ! Choice of boundary condition for sea ice variables 
-   INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim2_dta  !: = 0 use the initial state as bdy dta ; 
+#if ( defined key_lim2 || defined key_lim3 )
+   CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim       ! Choice of boundary condition for sea ice variables 
+   INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim_dta   !: = 0 use the initial state as bdy dta ; 
                                                               !: = 1 read it in a NetCDF file
 #endif
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4333)
@@ -197,5 +197,5 @@
 
 #if defined key_lim2
-            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
+            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_frld ) THEN
@@ -649,5 +649,5 @@
 #if defined key_lim2
             ! sea ice
-            IF( nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN
+            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN
 
                IF( dta%ll_frld ) THEN
@@ -724,4 +724,5 @@
                ENDIF
 
+            ENDIF
 #endif
             ! Recalculate field counts
@@ -850,5 +851,5 @@
 #if defined key_lim2
          IF (nn_ice_lim(ib_bdy) .gt. 0) THEN
-            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN
+            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN
                ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) )
                ALLOCATE( dta_bdy(ib_bdy)%hicif(nblen(1)) )
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 4333)
@@ -45,4 +45,5 @@
 
    REAL(wp) ::   epsi20 = 1.e-20_wp  ! module constants
+   REAL(wp) ::   epsi10 = 1.e-10_wp  ! min area allowed by ice model
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -91,13 +92,18 @@
       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index      !!
 
+      INTEGER  ::   jpbound            ! 0 = incoming ice
+                                       ! 1 = outgoing ice
       INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices
-      INTEGER  ::   ji, jj             ! local scalar
+      INTEGER  ::   ji, jj, ii, ij     ! local scalar
       REAL(wp) ::   zwgt, zwgt1        ! local scalar
-      REAL(wp) ::   zinda, ztmelts
+      REAL(wp) ::   zinda, ztmelts, zdh
+
+      REAL(wp), PARAMETER  ::   zsal = 6.3    ! arbitrary salinity    for incoming ice
+      REAL(wp), PARAMETER  ::   ztem = 270.0  ! arbitrary temperature for incoming ice
+      REAL(wp), PARAMETER  ::   zage = 30.0   ! arbitrary age         for incoming ice
       !!------------------------------------------------------------------------------
       !
       IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs')
       !
-      !IF( kt /= nit000 ) THEN
       jgrd = 1      ! Everything is at T-points here
       !
@@ -170,64 +176,92 @@
             ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth 
             ht_s(ji,jj,jl) = ( ht_s(ji,jj,jl) * zwgt1 + dta%ht_s(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
+
+            ! -----------------
+            ! Pathological case
+            ! -----------------
+            ! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to 
+            ! very large transformation from snow to ice (see limthd_dh.F90)
+
+            ! Then, a) 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 )
+            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
+            !zdh = MAX( 0._wp, ht_s(ji,jj,jl) * rhosn / rhoic )
+
+            ! recompute ht_i, ht_s
+            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 ) 
+
          ENDDO
          CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
          CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
          CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
-
+      ENDDO
+      ! retrieve at_i
+      at_i(:,:) = 0._wp
+      DO jl = 1, jpl
+         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+      END DO
+
+      DO jl = 1, jpl
          DO jb = 1, idx%nblen(jgrd)
             ji    = idx%nbi(jb,jgrd)
             jj    = idx%nbj(jb,jgrd)
-            ! clem : condition on ice thickness depends on the ice velocity
-            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values 
-            IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) THEN
-               ht_i(ji,jj,jl) = ht_i(ji+1,jj  ,jl)
-               a_i (ji,jj,jl) = a_i (ji+1,jj  ,jl)
-               ht_s(ji,jj,jl) = ht_s(ji+1,jj  ,jl)
-            ENDIF
-            IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) THEN
-               ht_i(ji,jj,jl) = ht_i(ji-1,jj  ,jl)
-               a_i (ji,jj,jl) = a_i (ji-1,jj  ,jl)
-               ht_s(ji,jj,jl) = ht_s(ji-1,jj  ,jl)
-            ENDIF
-            IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) THEN
-               ht_i(ji,jj,jl) = ht_i(ji  ,jj+1,jl)
-               a_i (ji,jj,jl) = a_i (ji  ,jj+1,jl)
-               ht_s(ji,jj,jl) = ht_s(ji  ,jj+1,jl)
-            ENDIF
-            IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) THEN
-               ht_i(ji,jj,jl) = ht_i(ji  ,jj-1,jl)
-               a_i (ji,jj,jl) = a_i (ji  ,jj-1,jl)
-               ht_s(ji,jj,jl) = ht_s(ji  ,jj-1,jl)
-            ENDIF
-
-            zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - a_i(ji,jj,jl) ) ) ! 0 if no ice
-            ! --------------------
+
+            ! condition on ice thickness depends on the ice velocity
+            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
+            jpbound = 0; ii = ji; ij = jj;
+
+            IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj
+            IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj
+            IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1
+            IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1
+
+            zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ii,ij) + 0.01 ) ) ! 0 if no ice
+
+            ! concentration and thickness
+            a_i (ji,jj,jl) = a_i (ii,ij,jl) * zinda
+            ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * zinda
+            ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * zinda
+
             ! Ice and snow volumes
-            ! --------------------
             v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
             v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl)
 
-            ! -------------
-            ! Ice salinity
-            !---------------
-            sm_i(ji,jj,jl)   = zinda * bulk_sal  + ( 1.0 - zinda ) * s_i_min
+            SELECT CASE( jpbound )
+
+            CASE( 0 ) ! velocity is inward
+
+               ! Ice salinity, age, temperature
+               sm_i(ji,jj,jl)   = zinda * zsal  + ( 1.0 - zinda ) * s_i_min
+               o_i(ji,jj,jl)    = zinda * zage  + ( 1.0 - zinda )
+               t_su(ji,jj,jl)   = zinda * ztem  + ( 1.0 - zinda ) * ztem
+               DO jk = 1, nlay_s
+                  t_s(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt
+               END DO
+               DO jk = 1, nlay_i
+                  t_i(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt 
+                  s_i(ji,jj,jk,jl) = zinda * zsal + ( 1.0 - zinda ) * s_i_min
+               END DO
+               
+            CASE( 1 ) ! velocity is outward
+ 
+               ! Ice salinity, age, temperature
+               sm_i(ji,jj,jl)   = zinda * sm_i(ii,ij,jl)  + ( 1.0 - zinda ) * s_i_min
+               o_i(ji,jj,jl)    = zinda * o_i(ii,ij,jl)   + ( 1.0 - zinda )
+               t_su(ji,jj,jl)   = zinda * t_su(ii,ij,jl)  + ( 1.0 - zinda ) * rtt
+               DO jk = 1, nlay_s
+                  t_s(ji,jj,jk,jl) = zinda * t_s(ii,ij,jk,jl) + ( 1.0 - zinda ) * rtt
+               END DO
+               DO jk = 1, nlay_i
+                  t_i(ji,jj,jk,jl) = zinda * t_i(ii,ij,jk,jl) + ( 1.0 - zinda ) * rtt
+                  s_i(ji,jj,jk,jl) = zinda * s_i(ii,ij,jk,jl) + ( 1.0 - zinda ) * s_i_min
+               END DO
+
+            END SELECT
+
+            ! contents
             smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
-            
-            !----------
-            ! Ice age
-            !----------
-            o_i(ji,jj,jl)    = zinda * 1.0   + ( 1.0 - zinda )
             oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
-            
-            !------------------------------
-            ! Sea ice surface temperature
-            !------------------------------
-            t_su(ji,jj,jl)   = zinda * 270.0 + ( 1.0 - zinda ) * t_bo(ji,jj)
-            
-            !------------------------------------
-            ! Snow temperature and heat content
-            !------------------------------------
             DO jk = 1, nlay_s
-               t_s(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt
                ! Snow energy of melting
                e_s(ji,jj,jk,jl) = zinda * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
@@ -235,33 +269,22 @@
                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
                ! Multiply by volume, so that heat content in 10^9 Joules
-               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
-                    v_s(ji,jj,jl)  / nlay_s
-            END DO !jk
-            
-            !-----------------------------------------------
-            ! Ice salinities, temperature and heat content 
-            !-----------------------------------------------
+               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s
+            END DO
             DO jk = 1, nlay_i
-               t_i(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt 
-               s_i(ji,jj,jk,jl) = zinda * bulk_sal + ( 1.0 - zinda ) * s_i_min
-               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
-               
+               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K                  
                ! heat content per unit volume
                e_i(ji,jj,jk,jl) = zinda * rhoic * &
-                    (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                    +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
-                    - rcp      * ( ztmelts - rtt ) )
-               
+                  (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                  +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
+                  - rcp      * ( ztmelts - rtt ) )
                ! Correct dimensions to avoid big values
                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
-               
                ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
-               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
-                    area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i
-            END DO ! jk
+               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i
+            END DO
+
 
          END DO !jb
  
-
          CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions
          CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
@@ -288,5 +311,4 @@
 #endif
       !      
-      !ENDIF !nit000/=0
       IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs')
       !
@@ -294,18 +316,16 @@
 
 
-   SUBROUTINE bdy_ice_lim_dyn( kn, pvar1, pvar2, pvar12 )
+   SUBROUTINE bdy_ice_lim_dyn( cd_type )
       !!------------------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_ice_lim_dyn  ***
       !!                    
       !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries.
-      !!      kn = 1: u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free
+      !!              u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free
       !!              if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities
-      !!      kn = 2: the stress tensor is set to 0 (i.e. pvar1, pvar2, pvar12)
       !!
       !! 2013-06 : C. Rousset
       !!------------------------------------------------------------------------------
       !!
-      INTEGER,  INTENT( in    )                           ::   kn     ! set up of ice vel (kn=1) or stress tensor (kn=2)
-      REAL(wp), INTENT( inout ), DIMENSION(:,:), OPTIONAL ::   pvar1, pvar2, pvar12  ! stress tensor components (zs1,zs2,zs12) 
+      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points
       INTEGER  ::   jb, jgrd   ! dummy loop indices
       INTEGER  ::   ji, jj             ! local scalar
@@ -317,88 +337,84 @@
       !
       DO ib_bdy=1, nb_bdy
-      !
-      SELECT CASE( nn_ice_lim(ib_bdy) )
-
-      CASE(jp_none)
-         CYCLE
-
-      CASE(jp_frs)
-
-         IF ( kn == 1 ) THEN ! set up of u_ice and v_ice
-           
-            jgrd = 2      ! u velocity
-            DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
-               ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
-               jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
-               zflag = idx_bdy(ib_bdy)%flagu(jb)
-
-               IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries
-                  ! one of the two zmsk is always 0 (because of zflag)
-                  zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
-                  zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
+         !
+         SELECT CASE( nn_ice_lim(ib_bdy) )
+
+         CASE('none')
+
+            CYCLE
+            
+         CASE('frs')
+            
+
+            SELECT CASE ( cd_type )
+
+            CASE ( 'U' )
+               
+               jgrd = 2      ! u velocity
+               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
+                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
+                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
+                  zflag = idx_bdy(ib_bdy)%flagu(jb)
                   
-                  ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
-                  u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
-                     &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
-                     &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
-               ELSE                             ! everywhere else
-                  u_ice(ji,jj) = u_oce(ji,jj)
-               ENDIF
-               ! mask ice velocities
-               zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice
-               u_ice(ji,jj) = zinda * u_ice(ji,jj)
-
-            ENDDO
+                  IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries
+                     ! one of the two zmsk is always 0 (because of zflag)
+                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
+                     
+                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
+                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
+                        &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
+                        &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
+                  ELSE                             ! everywhere else
+                     !u_ice(ji,jj) = u_oce(ji,jj)
+                     u_ice(ji,jj) = 0._wp
+                  ENDIF
+                  ! mask ice velocities
+                  zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice
+                  u_ice(ji,jj) = zinda * u_ice(ji,jj)
+                  
+               ENDDO
+
+               CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy )
+               
+            CASE ( 'V' )
+               
+               jgrd = 3      ! v velocity
+               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
+                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
+                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
+                  zflag = idx_bdy(ib_bdy)%flagv(jb)
+                  
+                  IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries
+                     ! one of the two zmsk is always 0 (because of zflag)
+                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
+                     
+                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
+                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
+                        &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
+                        &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
+                  ELSE                             ! everywhere else
+                     !v_ice(ji,jj) = v_oce(ji,jj)
+                     v_ice(ji,jj) = 0._wp
+                  ENDIF
+                  ! mask ice velocities
+                  zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice
+                  v_ice(ji,jj) = zinda * v_ice(ji,jj)
+                  
+               ENDDO
+               
+               CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy )
+               
+            END SELECT
             
-            jgrd = 3      ! v velocity
-            DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
-               ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
-               jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
-               zflag = idx_bdy(ib_bdy)%flagv(jb)
-              
-               IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries
-                  ! one of the two zmsk is always 0 (because of zflag)
-                  zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
-                  zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
-                  
-                  ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
-                  v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
-                     &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
-                     &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
-               ELSE                             ! everywhere else
-                  v_ice(ji,jj) = v_oce(ji,jj)
-               ENDIF
-               ! mask ice velocities
-               zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice
-               v_ice(ji,jj) = zinda * v_ice(ji,jj)
-
-            ENDDO
-            
-         CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy )
-         CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy )
- 
-         ELSEIF ( kn == 2 ) THEN ! set up of stress tensor (not sure it works)
+         CASE DEFAULT
+            CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' )
+         END SELECT
          
-            jgrd = 1      ! T grid
-            DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
-               ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
-               jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
-
-               pvar1 (ji,jj) = 0._wp
-               pvar2 (ji,jj) = 0._wp
-               pvar12(ji,jj) = 0._wp
-
-            ENDDO
-
-         ENDIF
-
-      CASE DEFAULT
-         CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' )
-      END SELECT
-
       ENDDO
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn')
- 
+      
     END SUBROUTINE bdy_ice_lim_dyn
 
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4333)
@@ -488,6 +488,8 @@
             DO igrd = 1, jpbgrd
                id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
-               nblendta(igrd,ib_bdy) = kdimsz(1)
-               jpbdtau = MAX(jpbdtau, kdimsz(1))
+               !clem nblendta(igrd,ib_bdy) = kdimsz(1)
+               !clem jpbdtau = MAX(jpbdtau, kdimsz(1))
+               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
+               jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))
             ENDDO
             CALL iom_close( inum )
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4333)
@@ -36,8 +36,8 @@
    LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
 
-   REAL(dp), SAVE                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
-   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ssh_ini              !
-   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
-   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcssh_loc_ini, scssh_loc_ini     !
+   REAL(wp), SAVE                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ssh_ini              !
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcssh_loc_ini, scssh_loc_ini     !
 
    !! * Substitutions
@@ -67,12 +67,12 @@
       !!
       INTEGER    ::   jk                          ! dummy loop indice
-      REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
-      REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
-      REAL(dp)   ::   z_hc        , z_sc          ! heat and salt content
-      REAL(dp)   ::   z_v1        , z_v2          ! volume
-      REAL(dp)   ::   zdeltat                     !    -     -
-      REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
-      REAL(dp)   ::   z_frc_trd_v                 !    -     -
-      REAL(dp), POINTER, DIMENSION(:,:)   ::   zsurf              !
+      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
+      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
+      REAL(wp)   ::   z_hc        , z_sc          ! heat and salt content
+      REAL(wp)   ::   z_v1        , z_v2          ! volume
+      REAL(wp)   ::   zdeltat                     !    -     -
+      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
+      REAL(wp)   ::   z_frc_trd_v                 !    -     -
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsurf              !
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('dia_hsb')      
@@ -104,20 +104,17 @@
       ! 2a -  Content variations !
       ! ------------------------ !
-      zdiff_v2 = 0._dp
-      zdiff_hc = 0._dp
-      zdiff_sc = 0._dp
+      zdiff_v2 = 0._wp
+      zdiff_hc = 0._wp
+      zdiff_sc = 0._wp
       ! volume variation (calculated with ssh)
       zdiff_v1 = glob_sum( zsurf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
       DO jk = 1, jpkm1
          ! volume variation (calculated with scale factors)
-         zdiff_v2 = zdiff_v2 &
-            &  + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
+         zdiff_v2 = zdiff_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
          ! heat content variation
-         zdiff_hc = zdiff_hc &
-            &  + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
+         zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
             &                           - hc_loc_ini(:,:,jk) ) )
          ! salt content variation
-         zdiff_sc = zdiff_sc &
-            &  + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
+         zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
             &                           - sc_loc_ini(:,:,jk) ) )
       ENDDO
@@ -140,19 +137,16 @@
       ! 2b -  Content           !
       ! ----------------------- !
-      z_v2 = 0._dp
-      z_hc = 0._dp
-      z_sc = 0._dp
+      z_v2 = 0._wp
+      z_hc = 0._wp
+      z_sc = 0._wp
       ! volume (calculated with ssh)
       z_v1 = glob_sum( zsurf(:,:) * sshn(:,:) )
       DO jk = 1, jpkm1
          ! volume (calculated with scale factors)
-         z_v2 = z_v2 &
-            &     + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
+         z_v2 = z_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
          ! heat content
-         z_hc = z_hc &
-            &     + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) )
+         z_hc = z_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) )
          ! salt content
-         z_sc = z_sc &
-            &     + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) )
+         z_sc = z_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) )
       ENDDO
       ! add ssh if not vvl
@@ -170,5 +164,5 @@
       CALL iom_put( 'bgtemper' , z_hc / z_v2 )                      ! Temperature (C) 
       CALL iom_put( 'bgsaline' , z_sc / z_v2 )                      ! Salinity (psu)
-      CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_dp ) ! Heat content variation (10^9 J)
+      CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J)
       CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9 )                 ! Salt content variation (psu*km3) 
       CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 )                    ! volume ssh (km3)  
@@ -176,5 +170,5 @@
       CALL iom_put( 'bgvoltot' , zdiff_v2 * 1.e-9 )                 ! volume total (km3) 
       CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 )                     ! vol - surface forcing (volume) 
-      CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_dp ) ! hc  - surface forcing (heat content) 
+      CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_wp ) ! hc  - surface forcing (heat content) 
       CALL iom_put( 'bgfrcsal' , frc_s * 1.e-9 )                     ! sc  - surface forcing (salt content) 
       !
@@ -224,4 +218,6 @@
       IF(lwp) THEN                   ! Control print
          WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
          WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
@@ -308,7 +304,7 @@
           hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
           scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
-          frc_v = 0._dp                                           
-          frc_t = 0._dp                                           
-          frc_s = 0._dp                                                  
+          frc_v = 0._wp                                           
+          frc_t = 0._wp                                           
+          frc_s = 0._wp                                                  
        ENDIF
 
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4333)
@@ -23,4 +23,5 @@
    USE eosbn2          ! equation of state            (eos bn2 routine)
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
+   USE domvvl          ! variable volume
    USE divcur          ! hor. divergence and curl      (div & cur routines)
    USE sbc_ice, ONLY : lk_lim3
@@ -210,4 +211,5 @@
          CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   )
          CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    )
+         IF( lk_vvl )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
       ELSE
          neuler = 0
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4333)
@@ -793,9 +793,9 @@
         Cd_n10(:,:) =   cdn_wave
       ELSE
-        Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
+        Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
       ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
-      Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b)
-      Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) !   L & Y eq. (6c), (6d)
+      Ce_n10  = 1.e-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b)
+      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) !   L & Y eq. (6c), (6d)
       !!
       !! Initializing transfert coefficients with their first guess neutral equivalents :
@@ -824,9 +824,9 @@
 
            !! Updating the neutral 10m transfer coefficients :
-           Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
+           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
            sqrt_Cd_n10 = sqrt(Cd_n10)
-           Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
+           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
            stab    = 0.5 + sign(0.5,zeta)
-           Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))           !  L & Y eq. (6c), (6d)
+           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab))           !  L & Y eq. (6c), (6d)
 
            !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) :
@@ -927,9 +927,9 @@
         Cd_n10(:,:) =   cdn_wave
       ELSE
-        Cd_n10  = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 
+        Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 
       ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
-      Ce_n10  = 1E-3*( 34.6 * sqrt_Cd_n10 )
-      Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1 - stab))
+      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 )
+      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
 
       !! Initializing transf. coeff. with their first guess neutral equivalents :
@@ -973,9 +973,9 @@
          ELSE
            !! Updating the neutral 10m transfer coefficients :
-           Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a)
+           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a)
            sqrt_Cd_n10 = sqrt(Cd_n10)
-           Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b)
+           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b)
            stab    = 0.5 + sign(0.5,zeta_u)
-           Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d)
+           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c-6d)
            !!
            !!
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4333)
@@ -58,4 +58,5 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
+   USE lib_fortran     ! 
 
 #if defined key_bdy 
@@ -167,6 +168,6 @@
          !
          IF( ln_nicep ) THEN      ! control print at a given point
-            jiindx = 6   ;   jjindx = 47
-            WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
+            jiindx = 177   ;   jjindx = 112
+            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
          ENDIF
       ENDIF
@@ -303,5 +304,6 @@
          fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp
          fhmec  (:,:) = 0._wp   ;   
-         fmmec  (:,:) = 0._wp     
+         fmmec  (:,:) = 0._wp
+         fmmflx (:,:) = 0._wp     
          focea2D(:,:) = 0._wp
          fsup2D (:,:) = 0._wp
@@ -328,13 +330,23 @@
                           CALL lim_rst_opn( kt )     ! Open Ice restart file
          !
-         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
-         !
-         IF( .NOT. lk_c1d ) THEN                     ! Ice dynamics & transport (except in 1D case)
+         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
+         ! ----------------------------------------------
+         ! ice dynamics and transport (except in 1D case)
+         ! ----------------------------------------------
+         IF( .NOT. lk_c1d ) THEN
                           CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
                           CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
                           CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
-         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
+         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
                           CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
                           CALL lim_var_agg( 1 ) 
+#if defined key_bdy
+                          ! bdy ice thermo 
+                          CALL lim_var_glo2eqv            ! equivalent variables
+                          CALL bdy_ice_lim( kt )
+                          CALL lim_itd_me_zapsmall
+                          CALL lim_var_agg(1)
+         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print
+#endif
                           CALL lim_update1
          ENDIF
@@ -349,7 +361,14 @@
                           old_oa_i(:,:,:)  = oa_i(:,:,:)
                           old_smv_i(:,:,:) = smv_i (:,:,:)
-         !                                           ! Ice thermodynamics 
+ 
+         ! ----------------------------------------------
+         ! ice thermodynamic
+         ! ----------------------------------------------
                           CALL lim_var_glo2eqv            ! equivalent variables
                           CALL lim_var_agg(1)             ! aggregate ice categories
+                          ! previous lead fraction and ice volume for flux calculations
+                          pfrld(:,:)   = 1._wp - at_i(:,:)
+                          phicif(:,:)  = vt_i(:,:)
+                          !
                           CALL lim_var_bv                 ! bulk brine volume (diag)
                           CALL lim_thd( kt )              ! Ice thermodynamics 
@@ -357,22 +376,16 @@
                           oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
                           CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin)
-         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
+         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  !
-         !
-         !                                           ! Global variables update
                           CALL lim_var_agg( 1 )           ! requested by limupdate
-                          CALL lim_update2                 ! Global variables update
-#if defined key_bdy
-                          CALL lim_var_glo2eqv            !
-                          CALL bdy_ice_lim( kt )          ! clem modif: bdy ice
-                          CALL bdy_ice_lim_dyn( 1 )       !?? repeated from limrhg.F90
-#endif
+                          CALL lim_update2                ! Global variables update
+
                           CALL lim_var_glo2eqv            ! equivalent variables (outputs)
                           CALL lim_var_agg(2)             ! aggregate ice thickness categories
-         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 2, ' - Final state - ' )   ! control print
+         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print
          !
                           CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes
          !
-         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
+         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
          !
          !                                           ! Diagnostics and outputs 
@@ -387,5 +400,5 @@
                           CALL lim_var_glo2eqv            ! ???
          !
-         IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash
+         IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash
          !
       ENDIF                                    ! End sea-ice time step only
@@ -414,5 +427,5 @@
 
 
-   SUBROUTINE lim_ctl
+   SUBROUTINE lim_ctl( kt )
       !!-----------------------------------------------------------------------
       !!                   ***  ROUTINE lim_ctl *** 
@@ -420,4 +433,5 @@
       !! ** Purpose :   Alerts in case of model crash
       !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt      ! ocean time step
       INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices
       INTEGER  ::   inb_altests       ! number of alert tests (max 20)
@@ -439,11 +453,11 @@
             DO ji = 1, jpi
                IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN
-                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
-                  WRITE(numout,*) ' at_i     ', at_i(ji,jj)
-                  WRITE(numout,*) ' Point - category', ji, jj, jl
-                  WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl)
-                  WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl)
-                  WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
-                  WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
+                  !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
+                  !WRITE(numout,*) ' at_i     ', at_i(ji,jj)
+                  !WRITE(numout,*) ' Point - category', ji, jj, jl
+                  !WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl)
+                  !WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl)
+                  !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
+                  !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                ENDIF
@@ -459,5 +473,5 @@
          DO ji = 1, jpi
             IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN
-               CALL lim_prt_state( ji, jj, 2, ' ALERTE 3 :   Very thick ice ' )
+               !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' )
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
             ENDIF
@@ -470,17 +484,17 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 0.5  .AND.  &
+            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  &
                &  at_i(ji,jj) > 0._wp   ) THEN
-               CALL lim_prt_state( ji, jj, 1, ' ALERTE 4 :   Very fast ice ' )
-               WRITE(numout,*) ' ice strength             : ', strength(ji,jj)
-               WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj) 
-               WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj)
-               WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj) 
-               WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj)
-               WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj)
-               WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj)
-               WRITE(numout,*) ' sst                      : ', sst_m(ji,jj)
-               WRITE(numout,*) ' sss                      : ', sss_m(ji,jj)
-               WRITE(numout,*) 
+               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' )
+               !WRITE(numout,*) ' ice strength             : ', strength(ji,jj)
+               !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj) 
+               !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj)
+               !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj) 
+               !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj)
+               !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj)
+               !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj)
+               !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj)
+               !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj)
+               !WRITE(numout,*) 
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
             ENDIF
@@ -494,13 +508,13 @@
          DO ji = 1, jpi
             IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN 
-               CALL lim_prt_state( ji, jj, 1, ' ALERTE 6 :   Ice on continents ' )
-               WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
-               WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
-               WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
-               WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj)
-               WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj)
-               WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1)
-               WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj)
-               WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj)
+               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' )
+               !WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
+               !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
+               !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
+               !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj)
+               !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj)
+               !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1)
+               !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj)
+               !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj)
                !
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
@@ -516,9 +530,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-!!gm  test twice sm_i ...  ????  bug?
-               IF( ( ( ABS( sm_i(ji,jj,jl) ) < 0.5 )   .OR.  &
-                     ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) ) .AND. &
-                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN
-!                 CALL lim_prt_state(ji,jj,1, ' ALERTE 7 :   Very fresh ice ' )
+               IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN
+!                 CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' )
 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
@@ -541,5 +552,5 @@
                       ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. &
                              ( a_i(ji,jj,jl) > 0._wp ) ) THEN
-                  CALL lim_prt_state( ji, jj, 1, ' ALERTE 9 :   Wrong ice age ')
+                  !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ')
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                ENDIF
@@ -553,14 +564,14 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN
-               CALL lim_prt_state( ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
-               DO jl = 1, jpl
-                  WRITE(numout,*) ' Category no: ', jl
-                  WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)   
-                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
-                  WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)   
-                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
-                  WRITE(numout,*) ' '
-               END DO
+            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth
+               !CALL lim_prt_state( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
+               !DO jl = 1, jpl
+                  !WRITE(numout,*) ' Category no: ', jl
+                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)   
+                  !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
+                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)   
+                  !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
+                  !WRITE(numout,*) ' '
+               !END DO
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
             ENDIF
@@ -573,26 +584,26 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF(   ABS( qns(ji,jj) ) > 1500._wp  .AND.  at_i(ji,jj) > 0._wp  )  THEN
+            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN
                !
-               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
-               WRITE(numout,*) ' ji, jj    : ', ji, jj
-               WRITE(numout,*) ' qns       : ', qns(ji,jj)
-               WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
-               WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
-               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj)
-               WRITE(numout,*) ' qldif     : ', qldif(ji,jj)
-               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice
-               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice
-               WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj)
-               WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj)
-               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice
-               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice
-               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
-               WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj) 
-               WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj) 
-               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
-               WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj) 
+               !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
+               !WRITE(numout,*) ' ji, jj    : ', ji, jj
+               !WRITE(numout,*) ' qns       : ', qns(ji,jj)
+               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
+               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
+               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj)
+               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj)
+               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice
+               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice
+               !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj)
+               !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj)
+               !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice
+               !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice
+               !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
+               !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj) 
+               !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj) 
+               !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
+               !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj) 
                !
-               CALL lim_prt_state( ji, jj, 2, '   ')
+               !CALL lim_prt_state( kt, ji, jj, 2, '   ')
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                !
@@ -611,12 +622,12 @@
                DO ji = 1, jpi
                   ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt
-                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-6   &
+                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   &
                      &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN
-                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
-                     WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
-                     WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
-                     WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
-                     WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl)
-                     WRITE(numout,*) ' ztmelts : ', ztmelts
+                     !WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
+                     !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
+                     !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
+                     !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
+                     !WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl)
+                     !WRITE(numout,*) ' ztmelts : ', ztmelts
                      inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                   ENDIF
@@ -626,16 +637,26 @@
       END DO
 
-      ialert_id = 1                                 ! reference number of this alert
-      cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert
-      WRITE(numout,*)
-      WRITE(numout,*) ' All alerts at the end of ice model '
-      DO ialert_id = 1, inb_altests
-         WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! '
-      END DO
-      !
+      ! sum of the alerts on all processors
+      IF( lk_mpp ) THEN
+         DO ialert_id = 1, inb_altests
+            CALL mpp_sum(inb_alp(ialert_id))
+         END DO
+      ENDIF
+
+      ! print alerts
+      IF( lwp ) THEN
+         ialert_id = 1                                 ! reference number of this alert
+         cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert
+         WRITE(numout,*) ' time step ',kt
+         WRITE(numout,*) ' All alerts at the end of ice model '
+         DO ialert_id = 1, inb_altests
+            WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! '
+         END DO
+      ENDIF
+     !
    END SUBROUTINE lim_ctl
  
    
-   SUBROUTINE lim_prt_state( ki, kj, kn, cd1 )
+   SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 )
       !!-----------------------------------------------------------------------
       !!                   ***  ROUTINE lim_prt_state *** 
@@ -651,172 +672,183 @@
       !!                n : number of the option
       !!-------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt      ! ocean time step
       INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices
       CHARACTER(len=*), INTENT(in) ::   cd1           !
       !!
-      INTEGER :: jl
+      INTEGER :: jl, ji, jj
       !!-------------------------------------------------------------------
 
-      WRITE(numout,*) cd1             ! print title
-
-      !----------------
-      !  Simple state
-      !----------------
-
-      IF ( kn == 1 .OR. kn == -1 ) THEN
-         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
-         WRITE(numout,*) ' Simple state '
-         WRITE(numout,*) ' masks s,u,v   : ', tms(ki,kj), tmu(ki,kj), tmv(ki,kj)
-         WRITE(numout,*) ' lat - long    : ', gphit(ki,kj), glamt(ki,kj)
-         WRITE(numout,*) ' Time step     : ', numit
-         WRITE(numout,*) ' - Ice drift   '
-         WRITE(numout,*) '   ~~~~~~~~~~~ '
-         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj)
-         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj)
-         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1)
-         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj)
-         WRITE(numout,*) ' strength      : ', strength(ki,kj)
-         WRITE(numout,*)
-         WRITE(numout,*) ' - Cell values '
-         WRITE(numout,*) '   ~~~~~~~~~~~ '
-         WRITE(numout,*) ' cell area     : ', area(ki,kj)
-         WRITE(numout,*) ' at_i          : ', at_i(ki,kj)       
-         WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)       
-         WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)       
-         DO jl = 1, jpl
-            WRITE(numout,*) ' - Category (', jl,')'
-            WRITE(numout,*) ' a_i           : ', a_i(ki,kj,jl)
-            WRITE(numout,*) ' ht_i          : ', ht_i(ki,kj,jl)
-            WRITE(numout,*) ' ht_s          : ', ht_s(ki,kj,jl)
-            WRITE(numout,*) ' v_i           : ', v_i(ki,kj,jl)
-            WRITE(numout,*) ' v_s           : ', v_s(ki,kj,jl)
-            WRITE(numout,*) ' e_s           : ', e_s(ki,kj,1,jl)/1.0e9
-            WRITE(numout,*) ' e_i           : ', e_i(ki,kj,1:nlay_i,jl)/1.0e9
-            WRITE(numout,*) ' t_su          : ', t_su(ki,kj,jl)
-            WRITE(numout,*) ' t_snow        : ', t_s(ki,kj,1,jl)
-            WRITE(numout,*) ' t_i           : ', t_i(ki,kj,1:nlay_i,jl)
-            WRITE(numout,*) ' sm_i          : ', sm_i(ki,kj,jl)
-            WRITE(numout,*) ' smv_i         : ', smv_i(ki,kj,jl)
-            WRITE(numout,*)
-            WRITE(numout,*) ' Pathological case : ', patho_case(ki,kj,jl)
-         END DO
-      ENDIF
-      IF( kn == -1 ) THEN
-         WRITE(numout,*) ' Mechanical Check ************** '
-         WRITE(numout,*) ' Check what means ice divergence '
-         WRITE(numout,*) ' Total ice concentration ', at_i (ki,kj)
-         WRITE(numout,*) ' Total lead fraction     ', ato_i(ki,kj)
-         WRITE(numout,*) ' Sum of both             ', ato_i(ki,kj) + at_i(ki,kj)
-         WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ki,kj) + at_i(ki,kj) - 1.00
-      ENDIF
-
-
-      !--------------------
-      !  Exhaustive state
-      !--------------------
-
-      IF ( kn .EQ. 2 ) THEN
-         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
-         WRITE(numout,*) ' Exhaustive state '
-         WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj)
-         WRITE(numout,*) ' Time step ', numit
-         WRITE(numout,*) 
-         WRITE(numout,*) ' - Cell values '
-         WRITE(numout,*) '   ~~~~~~~~~~~ '
-         WRITE(numout,*) ' cell area     : ', area(ki,kj)
-         WRITE(numout,*) ' at_i          : ', at_i(ki,kj)       
-         WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)       
-         WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)       
-         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj)
-         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj)
-         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1)
-         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj)
-         WRITE(numout,*) ' strength      : ', strength(ki,kj)
-         WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ki,kj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ki,kj)
-         WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ki,kj)  , ' old_v_ice     : ', old_v_ice(ki,kj)  
-         WRITE(numout,*)
-
-         DO jl = 1, jpl
-              WRITE(numout,*) ' - Category (',jl,')'
-              WRITE(numout,*) '   ~~~~~~~~         ' 
-              WRITE(numout,*) ' ht_i       : ', ht_i(ki,kj,jl)             , ' ht_s       : ', ht_s(ki,kj,jl)
-              WRITE(numout,*) ' t_i        : ', t_i(ki,kj,1:nlay_i,jl)
-              WRITE(numout,*) ' t_su       : ', t_su(ki,kj,jl)             , ' t_s        : ', t_s(ki,kj,1,jl)
-              WRITE(numout,*) ' sm_i       : ', sm_i(ki,kj,jl)             , ' o_i        : ', o_i(ki,kj,jl)
-              WRITE(numout,*) ' a_i        : ', a_i(ki,kj,jl)              , ' old_a_i    : ', old_a_i(ki,kj,jl)   
-              WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ki,kj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ki,kj,jl) 
-              WRITE(numout,*) ' v_i        : ', v_i(ki,kj,jl)              , ' old_v_i    : ', old_v_i(ki,kj,jl)   
-              WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ki,kj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ki,kj,jl) 
-              WRITE(numout,*) ' v_s        : ', v_s(ki,kj,jl)              , ' old_v_s    : ', old_v_s(ki,kj,jl)  
-              WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ki,kj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ki,kj,jl)
-              WRITE(numout,*) ' e_i1       : ', e_i(ki,kj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ki,kj,1,jl)/1.0e9 
-              WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ki,kj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ki,kj,1,jl)/1.0e9
-              WRITE(numout,*) ' e_i2       : ', e_i(ki,kj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ki,kj,2,jl)/1.0e9  
-              WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ki,kj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ki,kj,2,jl)/1.0e9
-              WRITE(numout,*) ' e_snow     : ', e_s(ki,kj,1,jl)            , ' old_e_snow : ', old_e_s(ki,kj,1,jl) 
-              WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ki,kj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ki,kj,1,jl)
-              WRITE(numout,*) ' smv_i      : ', smv_i(ki,kj,jl)            , ' old_smv_i  : ', old_smv_i(ki,kj,jl)   
-              WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ki,kj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ki,kj,jl) 
-              WRITE(numout,*) ' oa_i       : ', oa_i(ki,kj,jl)             , ' old_oa_i   : ', old_oa_i(ki,kj,jl)
-              WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ki,kj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ki,kj,jl)
-              WRITE(numout,*) ' Path. case : ', patho_case(ki,kj,jl)
-        END DO !jl
-
-        WRITE(numout,*)
-        WRITE(numout,*) ' - Heat / FW fluxes '
-        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
-!       WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ki,kj)
-!       WRITE(numout,*) ' sfx        : ', sfx      (ki,kj)
-!       WRITE(numout,*) ' sfx_res  : ', sfx_res(ki,kj)
-        WRITE(numout,*) ' fmmec      : ', fmmec    (ki,kj)
-        WRITE(numout,*) ' fhmec      : ', fhmec    (ki,kj)
-        WRITE(numout,*) ' fhbri      : ', fhbri    (ki,kj)
-        WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ki,kj)
-        WRITE(numout,*) 
-        WRITE(numout,*) ' sst        : ', sst_m(ki,kj)  
-        WRITE(numout,*) ' sss        : ', sss_m(ki,kj)  
-        WRITE(numout,*) 
-        WRITE(numout,*) ' - Stresses '
-        WRITE(numout,*) '   ~~~~~~~~ '
-        WRITE(numout,*) ' utau_ice   : ', utau_ice(ki,kj) 
-        WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ki,kj)
-        WRITE(numout,*) ' utau       : ', utau    (ki,kj) 
-        WRITE(numout,*) ' vtau       : ', vtau    (ki,kj)
-        WRITE(numout,*) ' oc. vel. u : ', u_oce   (ki,kj)
-        WRITE(numout,*) ' oc. vel. v : ', v_oce   (ki,kj)
-     ENDIF
-
-     !---------------------
-     ! Salt / heat fluxes
-     !---------------------
-
-     IF ( kn .EQ. 3 ) THEN
-        WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
-        WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
-        WRITE(numout,*) ' - Salt / Heat Fluxes '
-        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
-        WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj)
-        WRITE(numout,*) ' Time step ', numit
-        WRITE(numout,*)
-        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
-        WRITE(numout,*) ' qsr       : ', qsr(ki,kj)
-        WRITE(numout,*) ' qns       : ', qns(ki,kj)
-        WRITE(numout,*)
-        WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
-        WRITE(numout,*) ' emp       : ', emp    (ki,kj)
-        WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ki,kj)
-        WRITE(numout,*) ' sfx       : ', sfx    (ki,kj)
-        WRITE(numout,*) ' sfx_res   : ', sfx_res(ki,kj)
-        WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ki,kj)
-        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
-        WRITE(numout,*) ' fheat_res : ', fheat_res(ki,kj)
-        WRITE(numout,*)
-        WRITE(numout,*) ' - Momentum fluxes '
-        WRITE(numout,*) ' utau      : ', utau(ki,kj) 
-        WRITE(numout,*) ' vtau      : ', vtau(ki,kj)
-      ENDIF
-      WRITE(numout,*) ' '
-      !
+      DO ji = mi0(ki), mi1(ki)
+         DO jj = mj0(kj), mj1(kj)
+
+            WRITE(numout,*) ' time step ',kt,' ',cd1             ! print title
+
+            !----------------
+            !  Simple state
+            !----------------
+            
+            IF ( kn == 1 .OR. kn == -1 ) THEN
+               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj
+               WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
+               WRITE(numout,*) ' Simple state '
+               WRITE(numout,*) ' masks s,u,v   : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)
+               WRITE(numout,*) ' lat - long    : ', gphit(ji,jj), glamt(ji,jj)
+               WRITE(numout,*) ' Time step     : ', numit
+               WRITE(numout,*) ' - Ice drift   '
+               WRITE(numout,*) '   ~~~~~~~~~~~ '
+               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj)
+               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj)
+               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1)
+               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj)
+               WRITE(numout,*) ' strength      : ', strength(ji,jj)
+               WRITE(numout,*)
+               WRITE(numout,*) ' - Cell values '
+               WRITE(numout,*) '   ~~~~~~~~~~~ '
+               WRITE(numout,*) ' cell area     : ', area(ji,jj)
+               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
+               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
+               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)       
+               DO jl = 1, jpl
+                  WRITE(numout,*) ' - Category (', jl,')'
+                  WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl)
+                  WRITE(numout,*) ' ht_i          : ', ht_i(ji,jj,jl)
+                  WRITE(numout,*) ' ht_s          : ', ht_s(ji,jj,jl)
+                  WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl)
+                  WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl)
+                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9
+                  WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9
+                  WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl)
+                  WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl)
+                  WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl)
+                  WRITE(numout,*) ' sm_i          : ', sm_i(ji,jj,jl)
+                  WRITE(numout,*) ' smv_i         : ', smv_i(ji,jj,jl)
+                  WRITE(numout,*)
+               END DO
+            ENDIF
+            IF( kn == -1 ) THEN
+               WRITE(numout,*) ' Mechanical Check ************** '
+               WRITE(numout,*) ' Check what means ice divergence '
+               WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj)
+               WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj)
+               WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj)
+               WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00
+            ENDIF
+            
+
+            !--------------------
+            !  Exhaustive state
+            !--------------------
+            
+            IF ( kn .EQ. 2 ) THEN
+               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj
+               WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
+               WRITE(numout,*) ' Exhaustive state '
+               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
+               WRITE(numout,*) ' Time step ', numit
+               WRITE(numout,*) 
+               WRITE(numout,*) ' - Cell values '
+               WRITE(numout,*) '   ~~~~~~~~~~~ '
+               WRITE(numout,*) ' cell area     : ', area(ji,jj)
+               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
+               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
+               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)       
+               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj)
+               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj)
+               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1)
+               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj)
+               WRITE(numout,*) ' strength      : ', strength(ji,jj)
+               WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj)
+               WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ji,jj)  , ' old_v_ice     : ', old_v_ice(ji,jj)  
+               WRITE(numout,*)
+               
+               DO jl = 1, jpl
+                  WRITE(numout,*) ' - Category (',jl,')'
+                  WRITE(numout,*) '   ~~~~~~~~         ' 
+                  WRITE(numout,*) ' ht_i       : ', ht_i(ji,jj,jl)             , ' ht_s       : ', ht_s(ji,jj,jl)
+                  WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl)
+                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl)
+                  WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl)
+                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' old_a_i    : ', old_a_i(ji,jj,jl)   
+                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
+                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' old_v_i    : ', old_v_i(ji,jj,jl)   
+                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
+                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' old_v_s    : ', old_v_s(ji,jj,jl)  
+                  WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl)
+                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ji,jj,1,jl)/1.0e9 
+                  WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9
+                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ji,jj,2,jl)/1.0e9  
+                  WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9
+                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' old_e_snow : ', old_e_s(ji,jj,1,jl) 
+                  WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl)
+                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' old_smv_i  : ', old_smv_i(ji,jj,jl)   
+                  WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl) 
+                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' old_oa_i   : ', old_oa_i(ji,jj,jl)
+                  WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl)
+               END DO !jl
+               
+               WRITE(numout,*)
+               WRITE(numout,*) ' - Heat / FW fluxes '
+               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
+               WRITE(numout,*) ' emp        : ', emp      (ji,jj)
+               WRITE(numout,*) ' sfx        : ', sfx      (ji,jj)
+               WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj)
+               WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj)
+               WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj)
+               WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj)
+               WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj)
+               WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj)
+               WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj)
+               WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj)
+               WRITE(numout,*) 
+               WRITE(numout,*) ' sst        : ', sst_m(ji,jj)  
+               WRITE(numout,*) ' sss        : ', sss_m(ji,jj)  
+               WRITE(numout,*) 
+               WRITE(numout,*) ' - Stresses '
+               WRITE(numout,*) '   ~~~~~~~~ '
+               WRITE(numout,*) ' utau_ice   : ', utau_ice(ji,jj) 
+               WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ji,jj)
+               WRITE(numout,*) ' utau       : ', utau    (ji,jj) 
+               WRITE(numout,*) ' vtau       : ', vtau    (ji,jj)
+               WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj)
+               WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj)
+            ENDIF
+            
+            !---------------------
+            ! Salt / heat fluxes
+            !---------------------
+            
+            IF ( kn .EQ. 3 ) THEN
+               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj
+               WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
+               WRITE(numout,*) ' - Salt / Heat Fluxes '
+               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
+               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
+               WRITE(numout,*) ' Time step ', numit
+               WRITE(numout,*)
+               WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
+               WRITE(numout,*) ' qsr       : ', qsr(ji,jj)
+               WRITE(numout,*) ' qns       : ', qns(ji,jj)
+               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)
+               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice
+               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice
+               WRITE(numout,*)
+               WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
+               WRITE(numout,*) ' emp       : ', emp    (ji,jj)
+               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj)
+               WRITE(numout,*) ' sfx       : ', sfx    (ji,jj)
+               WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj)
+               WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj)
+               WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
+               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)
+               WRITE(numout,*)
+               WRITE(numout,*) ' - Momentum fluxes '
+               WRITE(numout,*) ' utau      : ', utau(ji,jj) 
+               WRITE(numout,*) ' vtau      : ', vtau(ji,jj)
+            ENDIF
+            WRITE(numout,*) ' '
+            !
+         END DO
+      END DO
+
    END SUBROUTINE lim_prt_state
 
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4333)
@@ -343,5 +343,5 @@
       !                                                           ! (update freshwater fluxes)
 !RBbug do not understand why see ticket 667
-      CALL lbc_lnk( emp, 'T', 1. )
+      !clem-bugsal CALL lbc_lnk( emp, 'T', 1. )
       !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 4332)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 4333)
@@ -33,5 +33,4 @@
    USE timing         ! Timing
    USE sbc_ice, ONLY : lk_lim3
-
 
    IMPLICIT NONE
