Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 6746)
@@ -80,5 +80,5 @@
       njeqm1 = njeq - 1 
 
-      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point
+      fcor(:,:) = 2. * omega * SIN( gphif(:,:) * rad )   !  coriolis factor at T-point
  
 !i    DO jj = 1, jpj
Index: anches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90	(revision 6745)
+++ 	(revision )
@@ -1,46 +1,0 @@
-MODULE dom_ice
-   !!======================================================================
-   !!                   ***  MODULE  dom_ice  ***
-   !! LIM-3 Sea Ice :   Domain  variables
-   !!======================================================================
-   !! History :  3.0  ! 2003-08  (M. Vancoppenolle)  LIM-3 original code
-   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation
-   !!----------------------------------------------------------------------
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC dom_ice_alloc   ! Routine called by nemogcm.F90
-
-   LOGICAL, PUBLIC ::   l_jeq = .TRUE.       !: Equator inside the domain flag
-
-   INTEGER, PUBLIC ::   njeq , njeqm1        !: j-index of the equator if it is inside the domain
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   fcor   !: coriolis coefficient
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   wght   !: weight of the 4 neighbours to compute averages
-
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   FUNCTION dom_ice_alloc()
-      !!-------------------------------------------------------------------
-      !!            *** Routine dom_ice_alloc ***
-      !!-------------------------------------------------------------------
-      INTEGER :: dom_ice_alloc
-      !!-------------------------------------------------------------------
-      !
-      ALLOCATE( fcor(jpi,jpj), wght(jpi,jpj,2,2), STAT = dom_ice_alloc )
-      !
-      IF( dom_ice_alloc /= 0 )   CALL ctl_warn( 'dom_ice_alloc: failed to allocate arrays.' )
-      !
-   END FUNCTION dom_ice_alloc
-
-   !!======================================================================
-END MODULE dom_ice
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6746)
@@ -203,5 +203,5 @@
    INTEGER , PUBLIC ::   nn_icestr        !: ice strength parameterization (0=Hibler79 1=Rothrock75)
    REAL(wp), PUBLIC ::   rn_pe_rdg        !: ridging work divided by pot. energy change in ridging, nn_icestr = 1
-   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength (N/M), Hibler JPO79
+   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79
    REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength
    LOGICAL , PUBLIC ::   ln_icestr_bvf    !: use brine volume to diminish ice strength
@@ -212,4 +212,7 @@
    INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
+   LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F) 
+   REAL(wp), PUBLIC ::   rn_gamma         !: fraction of ocean depth that ice must reach to initiate landfast ice
+   REAL(wp), PUBLIC ::   rn_icebfr        !: maximum bottom stress per unit area of contact (landfast ice) 
 
    !                                     !!** ice-diffusion namelist (namicehdf) **
@@ -296,5 +299,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field [s-1]
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sist        !: Average Sea-Ice Surface Temperature [Kelvin]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]     
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction
@@ -375,5 +377,4 @@
    !                                                                    !  this is an extensive variable that has to be transported
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age (days)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ov_i      !: Sea-Ice Age times volume per area (days.m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area (days)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i      !: brine volume
@@ -381,5 +382,4 @@
    !! Variables summed over all categories, or associated to all the ice in a single grid cell
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice !: components of the ice velocity (m/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tio_u, tio_v !: components of the ice-ocean stress (N/m2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s  !: ice and snow total volume per unit area (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i         !: ice total fractional area (ice concentration)
@@ -393,4 +393,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_s        !: mean snow thickness over all categories
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tau_icebfr   !: ice friction with bathy (landfast param activated)
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s      !: Snow temperatures [K]
@@ -466,12 +467,11 @@
       ! stay within Fortran's max-line length limit.
       ii = 1
-      ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ahiu      (jpi,jpj) , ahiv     (jpi,jpj) ,  &
-         &      hicol    (jpi,jpj) , strength (jpi,jpj) ,                                             &
-         &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                       &
-         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) , STAT=ierr(ii) )
-
-      ii = ii + 1
-      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                                           &
-         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        &
+      ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) ,                                             &
+         &      ahiu    (jpi,jpj) , ahiv     (jpi,jpj) , hicol    (jpi,jpj) ,                        &
+         &      strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  &
+         &      delta_i (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) , STAT=ierr(ii) )
+
+      ii = ii + 1
+      ALLOCATE( t_bo   (jpi,jpj) , frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,     &
          &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , wfx_lam(jpi,jpj) ,     &
          &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
@@ -493,10 +493,11 @@
          &      v_s    (jpi,jpj,jpl) , ht_s  (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,     &
          &      sm_i   (jpi,jpj,jpl) , smv_i (jpi,jpj,jpl) , o_i   (jpi,jpj,jpl) ,     &
-         &      ov_i   (jpi,jpj,jpl) , oa_i  (jpi,jpj,jpl) , bv_i  (jpi,jpj,jpl) ,  STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     &
+         &      oa_i   (jpi,jpj,jpl) , bv_i  (jpi,jpj,jpl) ,  STAT=ierr(ii) )
+      ii = ii + 1
+      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                       &
          &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     &
          &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     &
-         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) , om_i(jpi,jpj) , STAT=ierr(ii) )
+         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     &
+         &      om_i (jpi,jpj) , tau_icebfr(jpi,jpj)                              , STAT=ierr(ii) )
       ii = ii + 1
       ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
@@ -540,5 +541,5 @@
 
       ice_alloc = MAXVAL( ierr(:) )
-      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc_2: failed to allocate arrays.')
+      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc: failed to allocate arrays.')
       !
    END FUNCTION ice_alloc
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90	(revision 6746)
@@ -16,5 +16,4 @@
    !!----------------------------------------------------------------------
    USE dom_oce          ! ocean domain
-   USE dom_ice          ! LIM-3 domain
    USE ice              ! LIM-3 variables
    USE lbclnk           ! lateral boundary condition - MPP exchanges
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90	(revision 6746)
@@ -16,5 +16,4 @@
    USE dom_oce        ! ocean domain
    USE sbc_oce        ! ocean surface boundary condition
-   USE dom_ice        ! ice domain
    USE ice            ! ice variables
    !
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6746)
@@ -18,5 +18,4 @@
    USE phycst         ! physical constants
    USE ice            ! LIM-3 variables
-   USE dom_ice        ! LIM-3 domain
    USE dom_oce        ! ocean domain
    USE in_out_manager ! I/O manager
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 6746)
@@ -17,5 +17,4 @@
    USE ice             ! LIM-3: ice variables
    USE thd_ice         ! LIM-3: thermodynamical variables
-   USE dom_ice         ! LIM-3: ice domain
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice   fields
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 6746)
@@ -14,5 +14,4 @@
    !!----------------------------------------------------------------------
    USE ice             ! LIM-3: sea-ice variable
-   USE dom_ice         ! LIM-3: sea-ice domain
    USE dom_oce         ! ocean domain
    USE sbc_oce         ! surface boundary condition: ocean fields
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 6746)
@@ -19,5 +19,4 @@
    USE sbc_ice          ! Surface boundary condition: ice   fields
    USE ice              ! LIM-3 variables
-   USE dom_ice          ! LIM-3 domain
    USE limrhg           ! LIM-3 rheology
    USE lbclnk           ! lateral boundary conditions - MPP exchanges
@@ -58,9 +57,5 @@
       INTEGER, INTENT(in) ::   kt     ! number of iteration
       !!
-      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
-      INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
-      REAL(wp), POINTER, DIMENSION(:)   ::   zswitch        ! i-averaged indicator of sea-ice
-      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
-      !
+      INTEGER  :: jl, jk ! dummy loop indices
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
      !!---------------------------------------------------------------------
@@ -68,7 +63,5 @@
       IF( nn_timing == 1 )  CALL timing_start('limdyn')
 
-      CALL wrk_alloc( jpj, zswitch, zmsk )
-
-      CALL lim_var_agg(1)             ! aggregate ice categories
+      CALL lim_var_agg(1)                      ! aggregate ice categories
 
       IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
@@ -77,73 +70,23 @@
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       
+      ! ice velocities before rheology
       u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1)
       v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1)
       
+      ! Landfast ice parameterization: define max bottom friction
+      tau_icebfr(:,:) = 0._wp
+      IF( ln_landfast ) THEN
+         DO jl = 1, jpl
+            WHERE( ht_i(:,:,jl) > bathy(:,:) * rn_gamma )  tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
+         END DO
+      ENDIF
+      
       ! Rheology (ice dynamics)
-      ! ========
-      
-      !  Define the j-limits where ice rheology is computed
-      ! ---------------------------------------------------
-      
-      IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
-         i_j1 = 1
-         i_jpj = jpj
-         IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
-         CALL lim_rhg( i_j1, i_jpj )
-      ELSE                                 ! optimization of the computational area
-         !
-         DO jj = 1, jpj
-            zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
-            zmsk   (jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
-         END DO
-         
-         IF( l_jeq ) THEN                     ! local domain include both hemisphere
-            !                                 ! Rheology is computed in each hemisphere
-            !                                 ! only over the ice cover latitude strip
-            ! Northern hemisphere
-            i_j1  = njeq
-            i_jpj = jpj
-            DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
-               i_j1 = i_j1 + 1
-            END DO
-            i_j1 = MAX( 1, i_j1-2 )
-            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
-            CALL lim_rhg( i_j1, i_jpj )
-            !
-            ! Southern hemisphere
-            i_j1  =  1
-            i_jpj = njeq
-            DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
-               i_jpj = i_jpj - 1
-            END DO
-            i_jpj = MIN( jpj, i_jpj+1 )
-            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
-            !
-            CALL lim_rhg( i_j1, i_jpj )
-            !
-         ELSE                                 ! local domain extends over one hemisphere only
-            !                                 ! Rheology is computed only over the ice cover
-            !                                 ! latitude strip
-            i_j1  = 1
-            DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
-               i_j1 = i_j1 + 1
-            END DO
-            i_j1 = MAX( 1, i_j1-2 )
-            
-            i_jpj  = jpj
-            DO WHILE ( i_jpj >= 1  .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
-               i_jpj = i_jpj - 1
-            END DO
-            i_jpj = MIN( jpj, i_jpj+1)
-            !
-            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
-            !
-            CALL lim_rhg( i_j1, i_jpj )
-            !
-         ENDIF
-         !
-      ENDIF
+      ! ========     
+      CALL lim_rhg
       !
-      IF(ln_ctl) THEN   ! Control print
+      !
+      ! Control prints
+      IF(ln_ctl) THEN
          CALL prt_ctl_info(' ')
          CALL prt_ctl_info(' - Cell values : ')
@@ -173,10 +116,10 @@
             CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
             CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
-            DO ja = 1, nlay_i
+            DO jk = 1, nlay_i
                CALL prt_ctl_info(' ')
-               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
+               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
                CALL prt_ctl_info('   ~~~~~~~')
-               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
-               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
+               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_dyn  : t_i      : ')
+               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_dyn  : e_i      : ')
             END DO
          END DO
@@ -185,6 +128,4 @@
       ! conservation test
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      !
-      CALL wrk_dealloc( jpj, zswitch, zmsk )
       !
       IF( nn_timing == 1 )  CALL timing_stop('limdyn')
@@ -207,5 +148,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
-         &                nn_nevp, rn_relast
+         &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr
       !!-------------------------------------------------------------------
 
@@ -223,14 +164,17 @@
          WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr 
-         WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf
-         WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging   rn_pe_rdg     = ', rn_pe_rdg 
-         WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio
-         WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar
-         WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg
-         WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl
-         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc
-         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp
-         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast
+         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr 
+         WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf
+         WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging          rn_pe_rdg     = ', rn_pe_rdg 
+         WRITE(numout,*) '   drag coefficient for oceanic stress                         rn_cio        = ', rn_cio
+         WRITE(numout,*) '   first bulk-rheology parameter                               rn_pstar      = ', rn_pstar
+         WRITE(numout,*) '   second bulk-rhelogy parameter                               rn_crhg       = ', rn_crhg
+         WRITE(numout,*) '   creep limit                                                 rn_creepl     = ', rn_creepl
+         WRITE(numout,*) '   eccentricity of the elliptical yield curve                  rn_ecc        = ', rn_ecc
+         WRITE(numout,*) '   number of iterations for subcycling                         nn_nevp       = ', nn_nevp
+         WRITE(numout,*) '   ratio of elastic timescale over ice time step               rn_relast     = ', rn_relast
+         WRITE(numout,*) '   Landfast: param (T or F)                                    ln_landfast   = ', ln_landfast
+         WRITE(numout,*) '   Landfast: fraction of ocean depth that ice must reach       rn_gamma      = ', rn_gamma
+         WRITE(numout,*) '   Landfast: maximum bottom stress per unit area of contact    rn_icebfr     = ', rn_icebfr
       ENDIF
       !
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 6746)
@@ -23,5 +23,4 @@
    USE ice              ! sea-ice variables
    USE par_oce          ! ocean parameters
-   USE dom_ice          ! sea-ice domain
    USE limvar           ! lim_var_salprof
    USE in_out_manager   ! I/O manager
@@ -153,5 +152,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF( fcor(ji,jj) >= 0._wp ) THEN
+                  IF( ff(ji,jj) >= 0._wp ) THEN
                      zht_i_ini(ji,jj) = rn_hti_ini_n
                      zht_s_ini(ji,jj) = rn_hts_ini_n
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 6746)
@@ -18,5 +18,4 @@
    USE thd_ice          ! LIM thermodynamics
    USE ice              ! LIM variables
-   USE dom_ice          ! LIM domain
    USE limvar           ! LIM
    USE lbclnk           ! lateral boundary condition - MPP exchanges
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 6746)
@@ -18,5 +18,4 @@
    !!   lim_itd_shiftice :
    !!----------------------------------------------------------------------
-   USE dom_ice          ! LIM-3 domain
    USE par_oce          ! ocean parameters
    USE dom_oce          ! ocean domain
Index: anches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90	(revision 6745)
+++ 	(revision )
@@ -1,117 +1,0 @@
-MODULE limmsh
-   !!======================================================================
-   !!                     ***  MODULE  limmsh  ***
-   !! LIM ice model :   definition of the ice mesh parameters
-   !!======================================================================
-   !! History :  3.2  !  2008-01 (NEMO team)  LIM-3: adaptation from LIM-2
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_msh       : definition of the ice mesh
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constants
-   USE dom_oce        ! ocean domain
-   USE dom_ice        ! sea-ice domain
-   USE in_out_manager ! I/O manager
-   USE lbclnk         ! lateral boundary condition - MPP exchanges
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   lim_msh   ! routine called by sbcice_lim.F90
-
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_msh
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE lim_msh  ***
-      !!              
-      !! ** Purpose : Definition of the charact. of the numerical grid
-      !!       
-      !! ** Action  : - Initialisation of some variables
-      !!              - Definition of some constants linked with the grid
-      !!              - Definition of the metric coef. for the sea/ice
-      !! 
-      !! Reference  : Deleersnijder et al. Ocean Modelling 100, 7-10 
-      !!--------------------------------------------------------------------- 
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zusden   ! local scalar
-      !!---------------------------------------------------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_msh : LIM-3 sea-ice model, mesh initialization'
-         WRITE(numout,*) '~~~~~~~'
-      ENDIF
-
-      IF( jphgr_msh == 2 .OR. jphgr_msh == 3 .OR. jphgr_msh == 5 )   &
-          &      CALL ctl_stop(' Coriolis parameter in LIM not set for f- or beta-plane')
-
-      !                           !==  coriolis factor & Equator position ==!
-      njeq   = INT( jpj / 2 ) 
-      njeqm1 = njeq - 1 
-      !
-      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   ! coriolis factor
-      !
-      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere
-         l_jeq = .TRUE.
-         njeq  = 1
-         DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 )
-            njeq = njeq + 1
-         END DO
-         IF(lwp ) WRITE(numout,*) '          the equator is inside the domain at about njeq = ', njeq
-      ELSEIF( fcor(1,1) < 0.e0 ) THEN
-         l_jeq = .FALSE.
-         njeq = jpj
-         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the southern hemisphere: njeq = ', njeq
-      ELSE
-         l_jeq = .FALSE.
-         njeq = 2
-         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the northern hemisphere: njeq = ', njeq
-      ENDIF
-      !
-      njeqm1 = njeq - 1
-
-
-      !                           !==  metric coefficients for sea ice dynamic  ==!
-      wght(:,:,:,:) = 0.e0
-!!gm  Optimisation :  wght to be defined at F-point, not I-point  and change in limrhg
-      DO jj = 2, jpj
-         DO ji = 2, jpi
-            zusden = 1.e0 / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   &
-               &             * ( e2t(ji,jj) + e2t(ji  ,jj-1) ) )
-            wght(ji,jj,1,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
-            wght(ji,jj,1,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj-1)
-            wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj  )
-            wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1)
-         END DO
-      END DO
-      CALL lbc_lnk( wght(:,:,1,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point
-      CALL lbc_lnk( wght(:,:,1,2), 'I', 1. )      ! the value of wght at jpj is wrong
-      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
-      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
-!!gm end
-
-      !
-   END SUBROUTINE lim_msh
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option            Dummy Module         NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_msh           ! Dummy routine
-   END SUBROUTINE lim_msh
-#endif
-
-   !!======================================================================
-END MODULE limmsh
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6746)
@@ -9,10 +9,10 @@
    !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
    !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
-   !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
-   !!----------------------------------------------------------------------
-#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'               OR                     LIM-3 sea-ice model
-   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model
+   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
+   !!            3.6  !  2016-06  (C. Rousset) Rewriting and add the possibility to use mEVP from Bouillon 2013
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM-3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_rhg       : computes ice velocities
@@ -24,12 +24,6 @@
    USE sbc_oce        ! Surface boundary condition: ocean fields
    USE sbc_ice        ! Surface boundary condition: ice fields
-#if defined key_lim3
-   USE ice            ! LIM-3: ice variables
-   USE dom_ice        ! LIM-3: ice domain
-   USE limitd_me      ! LIM-3: 
-#else
-   USE ice_2          ! LIM-2: ice variables
-   USE dom_ice_2      ! LIM-2: ice domain
-#endif
+   USE ice            ! ice variables
+   USE limitd_me      ! ice strength
    USE lbclnk         ! Lateral Boundary Condition / MPP link
    USE lib_mpp        ! MPP library
@@ -38,8 +32,5 @@
    USE prtctl         ! Print control
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-#if defined key_agrif && defined key_lim2
-   USE agrif_lim2_interp
-#endif
-#if defined key_agrif && defined key_lim3
+#if defined key_agrif
    USE agrif_lim3_interp
 #endif
@@ -51,5 +42,5 @@
    PRIVATE
 
-   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
+   PUBLIC   lim_rhg        ! routine called by lim_dyn
 
    !! * Substitutions
@@ -62,5 +53,5 @@
 CONTAINS
 
-   SUBROUTINE lim_rhg( k_j1, k_jpj )
+   SUBROUTINE lim_rhg
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE lim_rhg  ***
@@ -109,133 +100,109 @@
       !!                 e.g. in the Canadian Archipelago
       !!
+      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
+      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
+      !!              but this solution appears very unstable (see Kimmritz et al 2016)
+      !!
       !! References : Hunke and Dukowicz, JPO97
       !!              Bouillon et al., Ocean Modelling 2009
+      !!              Bouillon et al., Ocean Modelling 2013
       !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
-      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
-      !!
-      INTEGER ::   ji, jj   ! dummy loop indices
-      INTEGER ::   jter     ! local integers
+      INTEGER ::   ji, jj       ! dummy loop indices
+      INTEGER ::   jter         ! local integers
       CHARACTER (len=50) ::   charout
-      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
-      REAL(wp) ::   za, zstms          ! local scalars
-      REAL(wp) ::   zc1, zc2, zc3      ! ice mass
-
-      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling
-      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity
-      REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars
-      REAL(wp) ::   zu_ice2, zv_ice1              !
-      REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre
-      REAL(wp) ::   zdst                          ! shear at the center of the grid point
-      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface
-      REAL(wp) ::   sigma1, sigma2                ! internal ice stress
-
-      REAL(wp) ::   zresm         ! Maximal error on ice velocity
-      REAL(wp) ::   zintb, zintn  ! dummy argument
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc)
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
-      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                           
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points
+
+      REAL(wp) ::   dtevp, z1_dtevp                            ! time step for subcycling
+      REAL(wp) ::   ecc2, z1_ecc2                              ! square of yield ellipse eccenticity
+      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2  ! alpha and beta from Bouillon 2009 and 2013
+      REAL(wp) ::   zm1, zm2, zm3                              ! ice/snow mass
+      REAL(wp) ::   delta, zp_delf, zds2
+      REAL(wp) ::   zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel  ! temporary scalars
+
+      REAL(wp) ::   zsig1, zsig2                               ! internal ice stress
+      REAL(wp) ::   zresm                                      ! Maximal error on ice velocity
+      REAL(wp) ::   zintb, zintn                               ! dummy argument
+
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                    ! temporary array for ice strength
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0          ! scale factors
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                   ! P/delta at T points
+     !
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za1   , za2               ! ice fraction on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2            ! ice/snow mass on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zcor1 , zcor2             ! coriolis parameter on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgu , zspgv             ! surface pressure gradient at U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1, u_oce2, v_ice1, u_ice2  ! ocean/ice u/v component on V/U points                           
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2               ! internal stresses
       
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2 
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope:
-                                                              !   ocean surface (ssh_m) if ice is not embedded
-                                                              !   ice top surface if ice is embedded   
-
-      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter
-      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt, zds                  ! tension and shear
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12            ! stress tensor components
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr     ! check convergence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                     ! array used for the calculation of ice surface slope:
+                                                                       !   ocean surface (ssh_m) if ice is not embedded
+                                                                       !   ice top surface if ice is embedded   
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitch1, zswitch2        ! dummy arrays
+
+      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp        ! tolerance parameter
+      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp        ! ice volume below which ice velocity equals ocean velocity
       !!-------------------------------------------------------------------
 
-      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
-      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
-      CALL wrk_alloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 )
-
-#if  defined key_lim2 && ! defined key_lim2_vp
-# if defined key_agrif
-      USE ice_2, vt_s => hsnm
-      USE ice_2, vt_i => hicm
-# else
-      vt_s => hsnm
-      vt_i => hicm
-# endif
-      at_i(:,:) = 1. - frld(:,:)
-#endif
-#if defined key_agrif && defined key_lim2 
-      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
-#endif
-#if defined key_agrif && defined key_lim3 
+      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
+      CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 )
+      CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 )
+      CALL wrk_alloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2 )
+
+#if defined key_agrif 
       CALL agrif_interp_lim3('U')   ! First interpolation of coarse values
-      CALL agrif_interp_lim3('V')   ! First interpolation of coarse values
-#endif
-      !
-      !------------------------------------------------------------------------------!
-      ! 1) Ice strength (zpresh)                                !
-      !------------------------------------------------------------------------------!
-      !
-      ! Put every vector to 0
-      delta_i(:,:) = 0._wp   ;
-      zpresh (:,:) = 0._wp   ;  
-      zpreshc(:,:) = 0._wp
-      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
-      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
-      shear_i(:,:) = 0._wp
-
-#if defined key_lim3
-      CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points
-#endif
-
-      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
-         DO ji = 1 , jpi
-#if defined key_lim3
-            zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj)
-#endif
-#if defined key_lim2
-            zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
-#endif
-            ! zmask = 1 where there is ice or on land
-            zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)
+      CALL agrif_interp_lim3('V')
+#endif
+      !
+      !------------------------------------------------------------------------------!
+      ! 0) define some variables
+      !------------------------------------------------------------------------------!
+      ! ecc2: square of yield ellipse eccenticrity
+      ecc2    = rn_ecc * rn_ecc
+      z1_ecc2 = 1._wp / ecc2
+
+      ! Time step for subcycling
+      dtevp    = rdt_ice / REAL( nn_nevp )
+      z1_dtevp = 1._wp / dtevp
+
+      ! alpha parameters (Bouillon 2009)
+      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+      zalph2 = zalph1 * z1_ecc2
+
+      ! alpha and beta parameters (Bouillon 2013)
+      !!zalph1 = 40.
+      !!zalph2 = 40.
+      !!zbeta  = 3000.
+      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
+
+      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+
+      !------------------------------------------------------------------------------!
+      ! 1) initialize arrays
+      !------------------------------------------------------------------------------!
+      ! Initialise stress tensor 
+      zs1 (:,:) = stress1_i (:,:) 
+      zs2 (:,:) = stress2_i (:,:)
+      zs12(:,:) = stress12_i(:,:)
+
+      ! Ice strength
+      CALL lim_itd_me_icestrength( nn_icestr )
+      zpresh(:,:) = tmask(:,:,1) *  strength(:,:)
+
+      ! scale factors
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
+            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
          END DO
       END DO
-
-      ! Ice strength on grid cell corners (zpreshc)
-      ! needed for calculation of shear stress 
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
-            zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   &
-               &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1)
-            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
-               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     &
-               &             ) / MAX( zstms, zepsi )
-         END DO
-      END DO
-      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
+            
       !
       !------------------------------------------------------------------------------!
       ! 2) Wind / ocean stress, mass terms, coriolis terms
       !------------------------------------------------------------------------------!
-      !
-      !  Wind stress, coriolis and mass terms on the sides of the squares        
-      !  zfrld1: lead fraction on U-points                                      
-      !  zfrld2: lead fraction on V-points                                     
-      !  zmass1: ice/snow mass on U-points                                    
-      !  zmass2: ice/snow mass on V-points                                   
-      !  zcorl1: Coriolis parameter on U-points                             
-      !  zcorl2: Coriolis parameter on V-points                            
-      !  (ztagnx,ztagny): wind stress on U/V points                       
-      !  v_oce1: ocean v component on u points                          
-      !  u_oce2: ocean u component on v points                         
 
       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
@@ -249,5 +216,5 @@
          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
          !
-         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
+         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
          !
       ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
@@ -255,52 +222,36 @@
       ENDIF
 
-      DO jj = k_j1+1, k_jpj-1
+      DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1
 
-            zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
-            zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
-            zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
-
-            zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj)
-            zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj)
-            zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  )
-            zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1)
-
-            ! Leads area.
-            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi )
-            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi )
-
-            ! Mass, coriolis coeff. and currents
-            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
-            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
-            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   &
-               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi )
-            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   &
-               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi )
-            !
-            ! Ocean has no slip boundary condition
-            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      &
-               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  &
-               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
-
-            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      &
-               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  &
-               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-
-            ! Wind stress at U,V-point
-            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
-            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
-
-            ! Computation of the velocity field taking into account the ice internal interaction.
-            ! Terms that are independent of the velocity field.
-
-            ! SB On utilise maintenant le gradient de la pente de l'ocean
-            ! include it later
-
-            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
-            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
-
-            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
-            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
+            ! ice fraction at U-V points
+            za1(ji,jj) = ( at_i(ji,jj) * e1t(ji+1,jj) + at_i(ji+1,jj) * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+            za2(ji,jj) = ( at_i(ji,jj) * e2t(ji,jj+1) + at_i(ji,jj+1) * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+
+            ! Ice/snow mass at U-V points
+            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
+            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
+            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
+            zmass1(ji,jj) = ( zm1 * e1t(ji+1,jj) + zm2 * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+            zmass2(ji,jj) = ( zm1 * e2t(ji,jj+1) + zm3 * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+
+            ! Ocean currents at U-V points
+            v_oce1(ji,jj)  = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
+               &                      + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+            
+            u_oce2(ji,jj)  = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
+               &                      + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+
+            ! Coriolis at U-V points (m*f)
+            zcor1(ji,jj) =   zmass1(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji  ,jj-1) )
+            zcor2(ji,jj) = - zmass2(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji-1,jj  ) )
+
+            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
+            zspgu(ji,jj) = - zmass1(ji,jj) * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            zspgv(ji,jj) = - zmass2(ji,jj) * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! switches
+            zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) )
+            zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) )
 
          END DO
@@ -312,264 +263,249 @@
       !------------------------------------------------------------------------------!
       !
-      ! Time step for subcycling
-      dtevp  = rdt_ice / nn_nevp
-#if defined key_lim3
-      dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )
-#else
-      dtotel = dtevp / ( 2._wp * telast )
-#endif
-      z1_dtotel = 1._wp / ( 1._wp + dtotel )
-      z1_dtevp  = 1._wp / dtevp
-      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
-      ecc2 = rn_ecc * rn_ecc
-      ecci = 1. / ecc2
-
-      !-Initialise stress tensor 
-      zs1 (:,:) = stress1_i (:,:) 
-      zs2 (:,:) = stress2_i (:,:)
-      zs12(:,:) = stress12_i(:,:)
-
       !                                               !----------------------!
       DO jter = 1 , nn_nevp                           !    loop over jter    !
          !                                            !----------------------!        
-         DO jj = k_j1, k_jpj-1
-            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
-            zv_ice(:,jj) = v_ice(:,jj)
-         END DO
-
-         DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask
-
-               !  
-               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
-               !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells
-               !- zds(:,:): shear on northeast corner of grid cells
-               !
-               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 
-               !                      there are many repeated calculations. 
-               !                      Speed could be improved by regrouping terms. For
-               !                      the moment, however, the stress is on clarity of coding to avoid
-               !                      bugs (Martin, for Miguel).
-               !
-               !- ALSO: arrays zdt, zds and delta could 
-               !  be removed in the future to minimise memory demand.
-               !
-               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
-               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
-               !              the corners is the same as in the B grid.
-               !
-               !
-               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+         ! Convergence test
+         IF(ln_ctl) THEN
+            DO jj = 1, jpjm1
+               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
+               zv_ice(:,jj) = v_ice(:,jj)
+            END DO
+         ENDIF
+
+         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               
+               ! divergence at T points
+               divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
                   &            ) * r1_e12t(ji,jj)
 
+               ! tension at T points
                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
                   &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
                   &         ) * r1_e12t(ji,jj)
 
-               !
+               ! shear at F points
                zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
                   &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-                  &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   &
-                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
-
-
-               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
-                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
-                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
-
-               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
-                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
-                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-            END DO
-         END DO
-
-         CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond.
+                  &         ) * r1_e12f(ji,jj)
+
+               ! u_ice at V point
+               u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+               
+               ! v_ice at U point
+               v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+
+            END DO
+         END DO
+         CALL lbc_lnk( zds, 'F', 1. )
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 ! no vector loop
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e12t(ji,jj)
+              
+               ! delta at T points
+               delta          = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zds2 ) * usecc2 )  
+               delta_i(ji,jj) = delta + rn_creepl
+
+               ! P/delta at T points
+               zp_delt(ji,jj) = zpresh(ji,jj) / delta_i(ji,jj)
+            END DO
+         END DO
+         CALL lbc_lnk( zp_delt, 'T', 1. )
+
+         ! --- Stress tensor --- !
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 ! no vector loop
+               
+               ! P/delta at F points
+               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
+               
+               ! stress tensor at T points
+               zs1(ji,jj) = ( zs1 (ji,jj) * zalph1 + zp_delt(ji,jj) * ( divu_i(ji,jj) - (delta_i(ji,jj)-rn_creepl) ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2 (ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt(ji,jj) * z1_ecc2  )                      ) * z1_alph2
+               ! F points
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf        * ( zds(ji,jj) * z1_ecc2  ) * 0.5_wp             ) * z1_alph2
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+ 
+         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1               
+               ! U points
+               zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e12u(ji,jj)
+
+               ! V points
+               zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e12v(ji,jj)
+            END DO
+         END DO
+         !
+         ! --- Computation of ice velocity --- !
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen smart and large nn_nevp
+         !  Bouillon et al. 2009 (eq 34-35) => stable
+         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
+            
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+
+                  zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
+
+                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj)
+                  zTauO   =   za2(ji,jj) * rhoco * zvel
+                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zCor    =   zcor2(ji,jj)  * u_ice2(ji,jj)
+                  zmt     =   zmass2(ji,jj) * z1_dtevp
+                  
+                  ! Bouillon 2009
+                  v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  &
+                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
+                  ! Bouillon 2013
+                  !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  &
+                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
+
+               END DO
+            END DO
+            CALL lbc_lnk( v_ice, 'V', -1. )
+            
+#if defined key_agrif
+            CALL agrif_interp_lim3( 'V' )
+#endif
+#if defined key_bdy
+            CALL bdy_ice_lim_dyn( 'V' )
+#endif         
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+
+                  ! v_ice at U point
+                  zv_ice1 = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+                     &               + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
+                  
+                  zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
+
+                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj)
+                  zTauO   =   za1(ji,jj) * rhoco * zvel
+                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zCor    =   zcor1(ji,jj)  * zv_ice1
+                  zmt     =   zmass1(ji,jj) * z1_dtevp
+                  
+                  ! Bouillon 2009
+                  u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  &
+                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
+                  ! Bouillon 2013
+                  !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!   &           + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  &
+                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( u_ice, 'U', -1. )
+            
+#if defined key_agrif
+            CALL agrif_interp_lim3( 'U' )
+#endif
+#if defined key_bdy
+            CALL bdy_ice_lim_dyn( 'U' )
+#endif         
+
+         ELSE ! odd iterations
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+
+                  zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) )
+
+                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj)
+                  zTauO   =   za1(ji,jj) * rhoco * zvel
+                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zCor    =   zcor1(ji,jj)  * v_ice1(ji,jj)
+                  zmt     =   zmass1(ji,jj) * z1_dtevp
+
+                  ! Bouillon 2009
+                  u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  &
+                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
+                  ! Bouillon 2013
+                  !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!   &           + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  &
+                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( u_ice, 'U', -1. )
+            
+#if defined key_agrif
+            CALL agrif_interp_lim3( 'U' )
+#endif
+#if defined key_bdy
+            CALL bdy_ice_lim_dyn( 'U' )
+#endif         
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+
+                  ! u_ice at V point
+                  zu_ice2 = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                     &               + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+
+                  zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) )
          
-         DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, fs_jpim1
-
-               !- Calculate Delta at centre of grid cells
-               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
-                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
-                  &            ) * r1_e12t(ji,jj)
-
-               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
-               delta_i(ji,jj) = delta + rn_creepl
-
-               !- Calculate Delta on corners
-               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
-                  &    ) * r1_e12f(ji,jj)
-
-               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
-                  &    ) * r1_e12f(ji,jj)
-
-               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl
-
-               !-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 * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   &
-                  &          ) * z1_dtotel
-               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   &
-                  &          ) * z1_dtotel
-               !-Calculate stress tensor component zs12 at corners
-               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  &
-                  &          ) * z1_dtotel 
-
-            END DO
-         END DO
-
-         CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
- 
-         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
-         DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, fs_jpim1
-               !- contribution of zs1, zs2 and zs12 to zf1
-               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
-                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          &
-                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  &
-                  &                ) * r1_e12u(ji,jj)
-               ! contribution of zs1, zs2 and zs12 to zf2
-               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
-                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          &
-                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  &
-                  &               )  * r1_e12v(ji,jj)
-            END DO
-         END DO
-         !
-         ! Computation of ice velocity
-         !
-         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
-         !
-         IF (MOD(jter,2).eq.0) THEN 
-
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = fs_2, fs_jpim1
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
-                  z0           = zmass1(ji,jj) * z1_dtevp
-
-                  ! SB modif because ocean has no slip boundary condition
-                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     &
-                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
-                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
-                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
-                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
-                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
+                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj)
+                  zTauO   =   za2(ji,jj) * rhoco * zvel
+                  ztauB   = - tau_icebfr(ji,jj) / zvel
+                  zCor    =   zcor2(ji,jj)  * zu_ice2
+                  zmt     =   zmass2(ji,jj) * z1_dtevp
+                  
+                  ! Bouillon 2009
+                  v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  &
+                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
+                  ! Bouillon 2013
+                  !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  &
+                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
+
                END DO
             END DO
-
-            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
-#if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
-#endif
-#if defined key_agrif && defined key_lim3
-            CALL agrif_interp_lim3( 'U' )
+            CALL lbc_lnk( v_ice, 'V', -1. )
+            
+#if defined key_agrif
+            CALL agrif_interp_lim3( 'V' )
 #endif
 #if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'U' )
-#endif         
-
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = fs_2, fs_jpim1
-
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
-                  z0           = zmass2(ji,jj) * z1_dtevp
-                  ! SB modif because ocean has no slip boundary condition
-                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
-                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
-                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
-                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) )
-                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
-               END DO
-            END DO
-
-            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
-#if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
-#endif
-#if defined key_agrif && defined key_lim3
-            CALL agrif_interp_lim3( 'V' )
-#endif
-#if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'V' )
-#endif         
-
-         ELSE 
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = fs_2, fs_jpim1
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
-                  z0           = zmass2(ji,jj) * z1_dtevp
-                  ! SB modif because ocean has no slip boundary condition
-                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
-                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
-                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)   
-
-                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &
-                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )
-                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
-               END DO
-            END DO
-
-            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
-#if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
-#endif
-#if defined key_agrif && defined key_lim3
-            CALL agrif_interp_lim3( 'V' )
-#endif
-#if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'V' )
-#endif         
-
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = fs_2, fs_jpim1
-                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
-                  z0           = zmass1(ji,jj) * z1_dtevp
-                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       &
-                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
-                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
-
-                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
-                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
-                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
-                  zcca         = z0 + za
-                  zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
-               END DO
-            END DO
-
-            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
-#if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
-#endif
-#if defined key_agrif && defined key_lim3
-            CALL agrif_interp_lim3( 'U' )
-#endif
-#if defined key_bdy
-         CALL bdy_ice_lim_dyn( 'U' )
+            CALL bdy_ice_lim_dyn( 'V' )
 #endif         
 
          ENDIF
          
+         ! Convergence test
          IF(ln_ctl) THEN
-            !---  Convergence test.
-            DO jj = k_j1+1 , k_jpj-1
+            DO jj = 2 , jpjm1
                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
             END DO
-            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
+            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
             IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
          ENDIF
-
+         !
          !                                                ! ==================== !
       END DO                                              !  end loop over jter  !
@@ -577,24 +513,19 @@
       !
       !------------------------------------------------------------------------------!
-      ! 4) Prevent ice velocities when the ice is thin
+      ! 4) Limit ice velocities when ice is thin
       !------------------------------------------------------------------------------!
       ! If the ice volume is below zvmin then ice velocity should equal the
       ! ocean velocity. This prevents high velocity when ice is thin
-      DO jj = k_j1+1, k_jpj-1
+       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1
-            IF ( vt_i(ji,jj) <= zvmin ) THEN
-               u_ice(ji,jj) = u_oce(ji,jj)
-               v_ice(ji,jj) = v_oce(ji,jj)
-            ENDIF
+            rswitch      = MAX( 0._wp, SIGN( 1._wp, vt_i(ji,jj) - zvmin ) )
+            u_ice(ji,jj) = ( rswitch * u_ice(ji,jj) + (1._wp - rswitch) * u_oce(ji,jj) ) * zswitch1(ji,jj)
+            v_ice(ji,jj) = ( rswitch * v_ice(ji,jj) + (1._wp - rswitch) * v_oce(ji,jj) ) * zswitch2(ji,jj)
          END DO
       END DO
 
-      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. )
-
-#if defined key_agrif && defined key_lim2
-      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' )
-      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' )
-#endif
-#if defined key_agrif && defined key_lim3
+      CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )
+
+#if defined key_agrif
       CALL agrif_interp_lim3( 'U' )
       CALL agrif_interp_lim3( 'V' )
@@ -605,73 +536,55 @@
 #endif         
 
-      DO jj = k_j1+1, k_jpj-1 
+      !------------------------------------------------------------------------------!
+      ! 5) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !------------------------------------------------------------------------------!
+
+      ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+      DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1
-            IF ( vt_i(ji,jj) <= zvmin ) THEN
-               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     &
-                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
-                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
-
-               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
-                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
-                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
-            ENDIF 
+            
+            ! divergence at T points
+            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+               &            ) * r1_e12t(ji,jj)
+            
+            ! tension at T points
+            zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+               &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               &         ) * r1_e12t(ji,jj)
+            
+            ! shear at F points
+            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e12f(ji,jj)
+            
          END DO
       END DO
-
-      CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. )
-
-      ! Recompute delta, shear and div, inputs for mechanical redistribution 
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask
-            !- divu_i(:,:), zdt(:,:): divergence and tension at centre 
-            !- zds(:,:): shear on northeast corner of grid cells
-            IF ( vt_i(ji,jj) <= zvmin ) THEN
-
-               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   &
-                  &            ) * r1_e12t(ji,jj)
-
-               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
-                  &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
-                  &         ) * r1_e12t(ji,jj)
-               !
-               ! SB modif because ocean has no slip boundary condition 
-               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
-                  &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     &
-                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
-
-               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    &
-                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj)
-
-               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
-               delta_i(ji,jj) = delta + rn_creepl
-            
-            ENDIF
+      CALL lbc_lnk_multi( divu_i, 'T', 1., zds, 'F', 1. )
+      
+
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! shear**2 at T points (doc eq. A16)
+            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
+               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
+               &   ) * 0.25_wp * r1_e12t(ji,jj)
+            
+            ! delta at T points
+            delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zds2 ) * usecc2 )  
+            delta_i(ji,jj) = delta + rn_creepl
+            
+            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zds2 )
          END DO
       END DO
-      !
-      !------------------------------------------------------------------------------!
-      ! 5) Store stress tensor and its invariants
-      !------------------------------------------------------------------------------!
-      ! * Invariants of the stress tensor are required for limitd_me
-      !   (accelerates convergence and improves stability)
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = fs_2, fs_jpim1
-            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &   
-               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 
-            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
-         END DO
-      END DO
-
-      ! Lateral boundary condition
-      CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. )
-
-      ! * Store the stress tensor for the next time step
+      CALL lbc_lnk_multi( delta_i, 'T', 1.,  shear_i, 'T', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
       stress1_i (:,:) = zs1 (:,:)
       stress2_i (:,:) = zs2 (:,:)
       stress12_i(:,:) = zs12(:,:)
-
-      !
+      !
+
       !------------------------------------------------------------------------------!
       ! 6) Control prints of residual and charge ellipse
@@ -695,9 +608,9 @@
             WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
             CALL prt_ctl_info(charout)
-            DO jj = k_j1+1, k_jpj-1
+            DO jj = 2, jpjm1
                DO ji = 2, jpim1
                   IF (zpresh(ji,jj) > 1.0) THEN
-                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
-                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
+                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) ) 
+                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) )
                      WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
                      CALL prt_ctl_info(charout)
@@ -710,8 +623,9 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
-      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
-      CALL wrk_dealloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 )
+      CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
+      CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 )
+      CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 )
+      CALL wrk_dealloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2 )
 
    END SUBROUTINE lim_rhg
@@ -722,6 +636,6 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
-      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
+   SUBROUTINE lim_rhg         ! Dummy routine
+      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?'
    END SUBROUTINE lim_rhg
 #endif
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 6746)
@@ -26,5 +26,4 @@
    USE sbc_ice        ! Surface boundary condition: ice fields
    USE thd_ice        ! LIM thermodynamic sea-ice variables
-   USE dom_ice        ! LIM sea-ice domain
    USE limthd_dif     ! LIM: thermodynamics, vertical diffusion
    USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 6746)
@@ -21,5 +21,4 @@
    USE sbc_ice        ! Surface boundary condition: ice fields
    USE thd_ice        ! LIM thermodynamics
-   USE dom_ice        ! LIM domain
    USE ice            ! LIM variables
    USE limtab         ! LIM 2D <==> 1D
@@ -71,8 +70,8 @@
       !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)      
       !!------------------------------------------------------------------------
-      INTEGER ::   ji,jj,jk,jl      ! dummy loop indices
-      INTEGER ::   nbpac            ! local integers 
-      INTEGER ::   ii, ij, iter     !   -       -
-      REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zde                         ! local scalars
+      INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices
+      INTEGER  ::   nbpac            ! local integers 
+      INTEGER  ::   ii, ij, iter     !   -       -
+      REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
@@ -154,6 +153,6 @@
 
       ! Default new ice thickness
-      WHERE( qlead(:,:) < 0._wp ) ; hicol = rn_hnewice
-      ELSEWHERE                   ; hicol = 0._wp
+      WHERE( qlead(:,:) < 0._wp ) ; hicol(:,:) = rn_hnewice
+      ELSEWHERE                   ; hicol(:,:) = 0._wp
       END WHERE
 
@@ -170,7 +169,7 @@
          zgamafr = 0.03
 
-         DO jj = 2, jpj
-            DO ji = 2, jpi
-               IF ( qlead(ji,jj) < 0._wp ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
                   !-------------
                   ! Wind stress
@@ -195,7 +194,6 @@
                   !-------------------
                   ! C-grid ice velocity
-                  rswitch = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
-                  zvgx    = rswitch * ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
-                  zvgy    = rswitch * ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
+                  zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
+                  zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
 
                   !-----------------------------------
@@ -203,6 +201,7 @@
                   !-----------------------------------
                   ! absolute relative velocity
-                  zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
-                     &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
+                  rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
+                  zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
+                     &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
                   zvrel(ji,jj) = SQRT( zvrel2 )
 
@@ -219,5 +218,5 @@
                      zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
 
-                     hicol(ji,jj) = hicol(ji,jj) - zf/zfp
+                     hicol(ji,jj) = hicol(ji,jj) - zf / MAX( zfp, epsi20 )
                      iter = iter + 1
                   END DO
@@ -228,6 +227,6 @@
          END DO 
          ! 
-         CALL lbc_lnk( zvrel(:,:), 'T', 1. )
-         CALL lbc_lnk( hicol(:,:), 'T', 1. )
+         CALL lbc_lnk( zvrel, 'T', 1. )
+         CALL lbc_lnk( hicol, 'T', 1. )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -240,5 +239,5 @@
       ! Select points for new ice formation
       !-------------------------------------
-      ! This occurs if open water energy budget is negative
+      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
       nbpac = 0
       npac(:) = 0
@@ -246,5 +245,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF ( qlead(ji,jj)  <  0._wp ) THEN
+            IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
                nbpac = nbpac + 1
                npac( nbpac ) = (jj - 1) * jpi + ji
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6746)
@@ -17,5 +17,4 @@
    USE dom_oce        ! ocean domain
    USE sbc_oce        ! ocean surface boundary condition
-   USE dom_ice        ! ice domain
    USE ice            ! ice variables
    USE limadv         ! ice advection
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 6746)
@@ -15,5 +15,4 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
-   USE dom_ice
    USE dom_oce
    USE phycst          ! physical constants
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 6746)
@@ -15,5 +15,4 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
-   USE dom_ice
    USE dom_oce
    USE phycst          ! physical constants
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6746)
@@ -41,5 +41,4 @@
    USE ice            ! ice variables
    USE thd_ice        ! ice variables (thermodynamics)
-   USE dom_ice        ! ice domain
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6746)
@@ -17,5 +17,4 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
-   USE dom_ice
    USE ice
    USE limvar
@@ -40,8 +39,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-#if defined key_dimgout
-# include "limwri_dimg.h90"
-#else
 
    SUBROUTINE lim_wri( kindic )
@@ -59,7 +54,7 @@
       INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::  z1_365
-      REAL(wp) ::  ztmp
+      REAL(wp) ::  z2da, z2db, ztmp
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace
       !!-------------------------------------------------------------------
 
@@ -67,5 +62,5 @@
 
       CALL wrk_alloc( jpi, jpj, jpl, zswi2 )
-      CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi )
+      CALL wrk_alloc( jpi, jpj     , z2d, zswi )
 
       !-----------------------------
@@ -95,18 +90,13 @@
          DO jj = 2 , jpjm1
             DO ji = 2 , jpim1
-               z2da(ji,jj)  = (  u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
-               z2db(ji,jj)  = (  v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
+               z2da  = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
+               z2db  = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
+               z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db )
            END DO
          END DO
-         CALL lbc_lnk( z2da, 'T', -1. )
-         CALL lbc_lnk( z2db, 'T', -1. )
-         CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component
-         CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component
-         DO jj = 1, jpj                                 
-            DO ji = 1, jpi
-               z2d(ji,jj)  = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) ) 
-            END DO
-         END DO
-         CALL iom_put( "icevel"       , z2d * zswi       )       ! ice velocity module
+         CALL lbc_lnk( z2d, 'T', 1. )
+         CALL iom_put( "uice_ipa"     , u_ice * zswi     )       ! ice velocity u component
+         CALL iom_put( "vice_ipa"     , v_ice * zswi     )       ! ice velocity v component
+         CALL iom_put( "icevel"       , z2d   * zswi     )       ! ice velocity module
       ENDIF
       !
@@ -130,5 +120,5 @@
       CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
 
-      CALL iom_put( "icestr"      , strength * 0.001 * zswi )    ! ice strength
+      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength
       CALL iom_put( "idive"       , divu_i * 1.0e8   * zswi )    ! divergence
       CALL iom_put( "ishear"      , shear_i * 1.0e8  * zswi )    ! shear
@@ -163,7 +153,16 @@
       CALL iom_put( "vfxlam"     , wfx_lam * ztmp       )        ! lateral melt 
       CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt 
+
+      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
+         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog
+         ELSEWHERE                                       ; z2d = 0._wp
+         END WHERE
+         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp )
+      ENDIF
+
+      ztmp = rday / rhosn
+      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
       CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt 
-      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow) 
-      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
+      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice) 
       
       CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total)
@@ -190,11 +189,4 @@
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
 
-
-      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
-         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog
-         ELSEWHERE                                       ; z2d = 0._wp
-         END WHERE
-         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp )
-      ENDIF
       
       !--------------------------------
@@ -223,10 +215,9 @@
       
       CALL wrk_dealloc( jpi, jpj, jpl, zswi2 )
-      CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db )
+      CALL wrk_dealloc( jpi, jpj     , z2d, zswi )
 
       IF( nn_timing == 1 )  CALL timing_stop('limwri')
       
    END SUBROUTINE lim_wri
-#endif
 
  
Index: anches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90	(revision 6745)
+++ 	(revision )
@@ -1,198 +1,0 @@
-SUBROUTINE lim_wri
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !!-------------------------------------------------------------------
-   !!  This routine computes the average of some variables and write it
-   !!  on the ouput files.
-   !!  ATTENTION cette routine n'est valable que si le pas de temps est
-   !!  egale a une fraction entiere de 1 jours.
-   !!  Diff 1-D 3-D : suppress common also included in etat
-   !!                 suppress cmoymo 11-18
-   !!  modif : 03/06/98
-   !!-------------------------------------------------------------------
-   USE  diawri, ONLY : dia_wri_dimg
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
-
-   REAL(wp),DIMENSION(1) ::   zdept
-
-   REAL(wp) ::   zsto, zsec, zjulian,zout, &
-   REAL(wp) ::   zindh,zinda,zindb, ztmu
-   REAL(wp), DIMENSION(jpi,jpj,jpnoumax) ::   zcmo
-   REAL(wp), DIMENSION(jpi,jpj) ::   zfield
-   INTEGER, SAVE ::   nmoyice   !: counter for averaging
-   INTEGER, SAVE ::   nwf       !: number of fields to write on disk
-   INTEGER, SAVE,DIMENSION (:), ALLOCATABLE  :: nsubindex   !: subindex to be saved
-   ! according to namelist
-
-   REAL(wp), SAVE, DIMENSION(jpi,jpj,jpnoumax) :: rcmoy
-#if ! defined key_diainstant
-   LOGICAL, PARAMETER :: ll_dia_inst=.false.      ! local logical variable 
-#else
-   LOGICAL, PARAMETER :: ll_dia_inst=.true.
-#endif
-   INTEGER ::  ji, jj, jf, ii   ! dummy loop indices and array index
-   INTEGER :: iyear, iday, imon ! 
-
-   CHARACTER(LEN=80) :: clname, cltext, clmode
-
-
-   INTEGER , SAVE ::   nice, nhorid, ndim, niter, ndepid
-   INTEGER , DIMENSION( jpij ) , SAVE ::   ndex51  
-   !!-------------------------------------------------------------------
-   IF ( numit == nstart ) THEN 
-
-      CALL lim_wri_init 
-
-      nwf = 0 
-      ii  = 0
-
-      IF (lwp ) THEN
-         WRITE(numout,*) 'lim_wri : Write ice outputs in dimg'
-         WRITE(numout,*) '~~~~~~~~'
-         WRITE(numout,*) '   According to namelist_ice, following fields saved:'
-         DO jf =1, noumef
-            IF (nc(jf) == 1 ) THEN
-               WRITE(numout,* ) '    -',titn(jf), nam(jf), uni(jf)
-            ENDIF
-         END DO
-      ENDIF
-
-      DO jf = 1, noumef
-         IF (nc(jf) == 1 ) nwf = nwf + 1
-      END DO
-
-      ALLOCATE( nsubindex (nwf) )
-
-      DO jf = 1, noumef
-         IF (nc(jf) == 1 ) THEN 
-            ii = ii +1 
-            nsubindex(ii) = jf
-         END IF
-      END DO
-
-      zsto     = rdt_ice
-      zout     = nwrite * rdt_ice / nn_fsbc
-      zsec     = 0.
-      niter    = 0
-      zdept(1) = 0.
-      nmoyice  = 0
-
-   ENDIF
-
-#if ! defined key_diainstant 
-   !-- calculs des valeurs instantanees
-
-   zcmo(:,:, 1:jpnoumax ) = 0.e0 
-   DO jj = 2 , jpjm1
-      DO ji = 2 , jpim1   ! NO vector opt.
-         zindh  = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,1) * (1.0 - frld(ji,jj) ) - 0.10 ) )
-         zinda  = MAX( 0._wp , SIGN( 1._wp , ( 1.0 - frld(ji,jj) ) - 0.10 ) )
-         zindb  = zindh * zinda
-         ztmu   = MAX( 0.5 * 1._wp , ( umask(ji,jj,1) + umask(ji+1,jj,1) + umask(ji,jj+1,1) + umask(ji+1,jj+1,1) ) ) 
-         zcmo(ji,jj,1)  = ht_s (ji,jj,1)
-         zcmo(ji,jj,2)  = ht_i (ji,jj,1)
-         zcmo(ji,jj,3)  = 0.
-         zcmo(ji,jj,4)  = frld  (ji,jj)
-         zcmo(ji,jj,5)  = sist  (ji,jj)
-         zcmo(ji,jj,6)  = fhtur  (ji,jj)
-         zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * umask(ji,jj,1) + u_ice(ji+1,jj  ) * umask(ji+1,jj,1)   &
-            + u_ice(ji,jj+1) * umask(ji,jj+1,1) + u_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) &
-            / ztmu 
-
-         zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * umask(ji,jj,1) + v_ice(ji+1,jj  ) * umask(ji+1,jj,1)   &
-            + v_ice(ji,jj+1) * umask(ji,jj+1,1) + v_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) &
-            / ztmu
-         zcmo(ji,jj,9)  = sst_m(ji,jj)
-         zcmo(ji,jj,10) = sss_m(ji,jj)
-
-         zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj)
-         zcmo(ji,jj,12) = qsr(ji,jj)
-         zcmo(ji,jj,13) = qns(ji,jj)
-         ! See thersf for the coefficient
-         zcmo(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in Kg/m2/day = mm/day
-         zcmo(ji,jj,15) = utau_ice(ji,jj)
-         zcmo(ji,jj,16) = vtau_ice(ji,jj)
-         zcmo(ji,jj,17) = qsr (ji,jj)
-         zcmo(ji,jj,18) = qns(ji,jj)
-         zcmo(ji,jj,19) = sprecip(ji,jj)
-      END DO
-   END DO
-   ! Cumulates values between outputs           
-   rcmoy(:,:,:)= rcmoy(:,:,:) + zcmo(:,:,:)
-   nmoyice = nmoyice + 1 
-   ! compute mean value if it is time to write on file
-   IF ( MOD(numit,nwrite) == 0 ) THEN
-      rcmoy(:,:,:) = rcmoy(:,:,:) / FLOAT(nmoyice)
-#else  
-      IF ( MOD(numit,nwrite) == 0 ) THEN 
-         !  case of instantaneaous output rcmoy(:,:, 1:jpnoumax ) = 0.e0
-         DO jj = 2 , jpjm1
-            DO ji = 2 , jpim1   ! NO vector opt.
-               zindh  = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,1) * (1.0 - frld(ji,jj) ) - 0.10 ) )
-               zinda  = MAX( 0._wp , SIGN( 1._wp , ( 1.0 - frld(ji,jj) ) - 0.10 ) )
-               zindb  = zindh * zinda
-               ztmu   = MAX( 0.5 * 1._wp , ( umask(ji,jj,1) + umask(ji+1,jj,1) + umask(ji,jj+1,1) + umask(ji+1,jj+1,1) ) )
-               rcmoy(ji,jj,1)  = ht_s (ji,jj,1)
-               rcmoy(ji,jj,2)  = ht_i (ji,jj,1)
-               rcmoy(ji,jj,3)  = 0.
-               rcmoy(ji,jj,4)  = frld  (ji,jj)
-               rcmoy(ji,jj,5)  = sist  (ji,jj)
-               rcmoy(ji,jj,6)  = fhtur  (ji,jj)
-               rcmoy(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * umask(ji,jj,1) + u_ice(ji+1,jj  ) * umask(ji+1,jj,1)   &
-                  + u_ice(ji,jj+1) * umask(ji,jj+1,1) + u_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) &
-                  / ztmu
-
-               rcmoy(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * umask(ji,jj,1) + v_ice(ji+1,jj  ) * umask(ji+1,jj,1)   &
-                  + v_ice(ji,jj+1) * umask(ji,jj+1,1) + v_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) &
-                  / ztmu
-               rcmoy(ji,jj,9)  = sst_m(ji,jj)
-               rcmoy(ji,jj,10) = sss_m(ji,jj)
-
-               rcmoy(ji,jj,11) = qns(ji,jj) + qsr(ji,jj)
-               rcmoy(ji,jj,12) = qsr(ji,jj)
-               rcmoy(ji,jj,13) = qns(ji,jj)
-               ! See thersf for the coefficient
-               rcmoy(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in mm/day
-               rcmoy(ji,jj,15) = utau_ice(ji,jj)
-               rcmoy(ji,jj,16) = vtau_ice(ji,jj)
-               rcmoy(ji,jj,17) = qsr(ji,jj)
-               rcmoy(ji,jj,18) = qns(ji,jj)
-               rcmoy(ji,jj,19) = sprecip(ji,jj)
-            END DO
-         END DO
-#endif
-
-         !
-         niter = niter + 1
-         DO jf = 1 , noumef
-            zfield(:,:) = (rcmoy(:,:,jf) * cmulti(jf) + cadd(jf)) * tmask(:,:,1)
-
-            IF ( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN 
-               CALL lbc_lnk( zfield, 'T', -1. )
-            ELSE 
-               CALL lbc_lnk( zfield, 'T',  1. )
-            ENDIF
-            rcmoy(:,:,jf) = zfield(:,:)
-         END DO
-
-         IF (ll_dia_inst) THEN
-            clmode='instantaneous'
-         ELSE
-            WRITE(clmode,'(f5.1,a)' ) nwrite*rdt/86400.,' days average'
-         END IF
-         iyear = ndastp/10000
-         imon = (ndastp-iyear*10000)/100
-         iday = ndastp - imon*100 - iyear*10000
-         WRITE(clname,9000) TRIM(cexper),'ICEMOD',iyear,imon,iday
-         cltext=TRIM(cexper)//' ice modele output'//TRIM(clmode)
-         CALL dia_wri_dimg (clname, cltext, rcmoy, nwf , 'I', nsubindex)
-9000     FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")
-
-         rcmoy(:,:,:) = 0.0
-         nmoyice = 0 
-      END IF     !  MOD(numit, nwrite == 0 ) !
-
-   END SUBROUTINE lim_wri
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 6746)
@@ -23,5 +23,4 @@
    USE sbc_oce
    USE ice
-   USE dom_ice
    USE agrif_ice
 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90	(revision 6746)
@@ -26,5 +26,4 @@
    USE agrif_oce
    USE ice
-   USE dom_ice
    USE agrif_ice 
 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 6746)
@@ -27,5 +27,4 @@
 #elif defined key_lim3
    USE ice             ! LIM_3 ice variables
-   USE dom_ice         ! sea-ice domain
    USE limvar
 #endif 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 6746)
@@ -211,21 +211,4 @@
       REAL(wp) ::   zztmp  
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
-      ! reading initial file
-      LOGICAL  ::   ln_tsd_init      !: T & S data flag
-      LOGICAL  ::   ln_tsd_tradmp    !: internal damping toward input data flag
-      CHARACTER(len=100)            ::   cn_dir
-      TYPE(FLD_N)                   ::  sn_tem,sn_sal
-      INTEGER  ::   ios=0
-
-      NAMELIST/namtsd/ ln_tsd_init,ln_tsd_tradmp,cn_dir,sn_tem,sn_sal
-      !
-
-      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :
-      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in reference namelist for dia_ar5', lwp )
-      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
-      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in configuration namelist for dia_ar5', lwp )
-      IF(lwm) WRITE ( numond, namtsd )
       !
       !!----------------------------------------------------------------------
@@ -233,5 +216,5 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
       !
-      CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_alloc( jpi, jpj, jpk, 2, zsaldta )
       !                                      ! allocate dia_ar5 arrays
       IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
@@ -249,8 +232,9 @@
       IF( lk_mpp )   CALL mpp_sum( vol0 )
 
-      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_sal%clname), inum )
-      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,1), 1  )
-      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 )
+      CALL iom_open ( 'sali_ref_clim_monthly', inum )
+      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  )
+      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
       CALL iom_close( inum )
+
       sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )        
       sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
@@ -267,5 +251,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_dealloc( jpi, jpj, jpk, 2, zsaldta )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 6746)
@@ -157,4 +157,14 @@
          END DO
       ENDIF
+
+      ! ORCA R1: Take the minimum between aeiw  and aeiv0
+      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 )
+            END DO
+         END DO
+      ENDIF
+
       CALL lbc_lnk( aeiw, 'W', 1. )      ! lateral boundary condition on aeiw 
 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6746)
@@ -1335,39 +1335,48 @@
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
       !!
-      !! ** Purpose :   provide the heat and freshwater fluxes of the 
-      !!              ocean-ice system.
+      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system
       !!
       !! ** Method  :   transform the fields received from the atmosphere into
       !!             surface heat and fresh water boundary condition for the 
       !!             ice-ocean system. The following fields are provided:
-      !!              * total non solar, solar and freshwater fluxes (qns_tot, 
+      !!               * total non solar, solar and freshwater fluxes (qns_tot, 
       !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
       !!             NB: emp_tot include runoffs and calving.
-      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
+      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
       !!             emp_ice = sublimation - solid precipitation as liquid
       !!             precipitation are re-routed directly to the ocean and 
-      !!             runoffs and calving directly enter the ocean.
-      !!              * solid precipitation (sprecip), used to add to qns_tot 
+      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90)
+      !!               * solid precipitation (sprecip), used to add to qns_tot 
       !!             the heat lost associated to melting solid precipitation
       !!             over the ocean fraction.
-      !!       ===>> CAUTION here this changes the net heat flux received from
-      !!             the atmosphere
-      !!
-      !!                  - the fluxes have been separated from the stress as
-      !!                 (a) they are updated at each ice time step compare to
-      !!                 an update at each coupled time step for the stress, and
-      !!                 (b) the conservative computation of the fluxes over the
-      !!                 sea-ice area requires the knowledge of the ice fraction
-      !!                 after the ice advection and before the ice thermodynamics,
-      !!                 so that the stress is updated before the ice dynamics
-      !!                 while the fluxes are updated after it.
+      !!               * heat content of rain, snow and evap can also be provided,
+      !!             otherwise heat flux associated with these mass flux are
+      !!             guessed (qemp_oce, qemp_ice)
+      !!
+      !!             - the fluxes have been separated from the stress as
+      !!               (a) they are updated at each ice time step compare to
+      !!               an update at each coupled time step for the stress, and
+      !!               (b) the conservative computation of the fluxes over the
+      !!               sea-ice area requires the knowledge of the ice fraction
+      !!               after the ice advection and before the ice thermodynamics,
+      !!               so that the stress is updated before the ice dynamics
+      !!               while the fluxes are updated after it.
+      !!
+      !! ** Details
+      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided
+      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
+      !!
+      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided
+      !!
+      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce)
+      !!                                                                      river runoff (rnf) is provided but not included here
       !!
       !! ** Action  :   update at each nf_ice time step:
       !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
       !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
-      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
-      !!                   emp_ice            ice sublimation - solid precipitation over the ice
-      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
-      !!                   sprecip             solid precipitation over the ocean  
+      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving)
+      !!                   emp_ice           ice sublimation - solid precipitation over the ice
+      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice
+      !!                   sprecip           solid precipitation over the ocean  
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
@@ -1379,5 +1388,5 @@
       INTEGER ::   jl         ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
@@ -1387,5 +1396,5 @@
       !
       CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
-      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
@@ -1396,29 +1405,27 @@
       !
       !                                                      ! ========================= !
-      !                                                      !    freshwater budget      !   (emp)
+      !                                                      !    freshwater budget      !   (emp_tot)
       !                                                      ! ========================= !
       !
-      !                                                           ! total Precipitation - total Evaporation (emp_tot)
-      !                                                           ! solid precipitation - sublimation       (emp_ice)
-      !                                                           ! solid Precipitation                     (sprecip)
-      !                                                           ! liquid + solid Precipitation            (tprecip)
+      !                                                           ! solid Precipitation                                (sprecip)
+      !                                                           ! liquid + solid Precipitation                       (tprecip)
+      !                                                           ! total Evaporation - total Precipitation            (emp_tot)
+      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice)
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
-      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
-         zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
-         ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
-         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
-         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
-            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation 
+      CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
+         zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
+         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
+         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
+         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:)
+               CALL iom_put( 'rain'         ,   frcv(jpr_rain)%z3(:,:,1)                                                         )  ! liquid precipitation 
          IF( iom_use('hflx_rain_cea') )   &
-            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip. 
-         IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
-            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
+            &  CALL iom_put( 'hflx_rain_cea',   frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:)                                            )  ! heat flux from liq. precip. 
          IF( iom_use('evap_ao_cea'  ) )   &
-            CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average)
+            &  CALL iom_put( 'evap_ao_cea'  ,   frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)                )  ! ice-free oce evap (cell average)
          IF( iom_use('hflx_evap_cea') )   &
-            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average)
-      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
+            &  CALL iom_put( 'hflx_evap_cea', ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) * zcptn(:,:) )  ! heat flux from from evap (cell average)
+      CASE( 'oce and ice' )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
-         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
+         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:)
          zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
@@ -1426,9 +1433,15 @@
 
 #if defined key_lim3
-      ! zsnw = snow percentage over ice after wind blowing
-      zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( p_frld, zsnw )
+      ! zsnw = snow fraction over ice after wind blowing
+      zsnw(:,:) = 0._wp  ;  CALL lim_thd_snwblow( p_frld, zsnw )
       
-      ! --- evaporation (kg/m2/s) --- !
+      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
+      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( zicefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
+      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
+
+      ! --- evaporation over ocean (used later for qemp) --- !
+      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
+
+      ! --- evaporation over ice (kg/m2/s) --- !
       zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
       ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
@@ -1436,14 +1449,11 @@
       zdevap_ice(:,:) = 0._wp
       
-      ! --- evaporation minus precipitation corrected for the effect of wind blowing on snow --- !
-      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) - zsprecip * (1._wp - zsnw)
-      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip * (1._wp - zsnw)          
-
-      ! Sublimation over sea-ice (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:) )
-      ! runoffs and calving (put in emp_tot)
+      ! --- runoffs (included in emp later on) --- !
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+
+      ! --- calving (put in emp_tot and emp_oce) --- !
       IF( srcv(jpr_cal)%laction ) THEN 
          zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
+         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
          CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
       ENDIF
@@ -1471,10 +1481,9 @@
       ENDIF
 
-                                     CALL iom_put( 'snowpre'    , sprecip                         )  ! Snow
-      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw ) )  ! Snow over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw   )  ! Snow over sea-ice         (cell average)    
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:)         )  ! Sublimation over sea-ice (cell average)
+                                     CALL iom_put( 'snowpre'    , sprecip(:,:)                         )  ! Snow
+      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw(:,:) ) )  ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw(:,:)   )  ! Snow over sea-ice         (cell average)
 #else
-      ! Sublimation over sea-ice (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )
       ! runoffs and calving (put in emp_tot)
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
@@ -1496,9 +1505,8 @@
       ENDIF
 
-         CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
-      IF( iom_use('snow_ao_cea') )   &
-         CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )   &
-         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )  ! Sublimation over sea-ice (cell average)
+                                    CALL iom_put( 'snowpre'    , sprecip(:,:)               )   ! Snow
+      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:) )   ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:) )   ! Snow over sea-ice         (cell average)
 #endif
 
@@ -1506,18 +1514,17 @@
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
       !                                                      ! ========================= !
-      CASE( 'oce only' )                                     ! the required field is directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
-      CASE( 'conservative' )                                      ! the required fields are directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
+      CASE( 'oce only' )         ! the required field is directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'conservative' )     ! the required fields are directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
-            ! Set all category values equal for the moment
             DO jl=1,jpl
-               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
             ENDDO
          ENDIF
-      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
-         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
+         zqns_tot(:,:) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             DO jl=1,jpl
@@ -1526,5 +1533,5 @@
             ENDDO
          ELSE
-            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+            qns_tot(:,:) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl=1,jpl
                zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
@@ -1532,10 +1539,10 @@
             ENDDO
          ENDIF
-      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
+      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
 ! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
          zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
             &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
-            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
+            &                                           + pist(:,:,1) * zicefr(:,:) ) )
       END SELECT
 !!gm
@@ -1547,18 +1554,11 @@
 !! similar job should be done for snow and precipitation temperature
       !                                     
-      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting 
-         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting 
-         zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:)
-         IF( iom_use('hflx_cal_cea') )   &
-            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
-      ENDIF
-
-      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus
-      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
+      IF( srcv(jpr_cal)%laction ) THEN   ! Iceberg melting 
+         zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! add the latent heat of iceberg melting
+                                                                         ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
+         IF( iom_use('hflx_cal_cea') )   CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus )   ! heat flux from calving
+      ENDIF
 
 #if defined key_lim3      
-      ! --- evaporation --- !
-      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
-
       ! --- non solar flux over ocean --- !
       !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
@@ -1567,5 +1567,5 @@
 
       ! --- heat flux associated with emp (W/m2) --- !
-      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &       ! evap
+      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn(:,:)   &       ! evap
          &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &       ! liquid precip
          &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus )  ! solid precip over ocean + snow melting
@@ -1606,4 +1606,12 @@
          qemp_ice (:,:  ) = zqemp_ice (:,:  )
       ENDIF
+
+      !! clem: we should output qemp_oce and qemp_ice (at least)
+      IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average)
+      !! these diags are not outputed yet
+!!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average)
+!!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average)
+!!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average)
+
 #else
       ! clem: this formulation is certainly wrong... but better than it was...
@@ -1611,5 +1619,5 @@
          &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
          &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
-         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
+         &             - zemp_ice(:,:) ) * zcptn(:,:) 
 
      IF( ln_mixcpl ) THEN
@@ -1731,5 +1739,5 @@
 
       CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
-      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
       CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6746)
@@ -24,5 +24,4 @@
    USE ice             ! LIM-3: ice variables
    USE thd_ice         ! LIM-3: thermodynamical variables
-   USE dom_ice         ! LIM-3: ice domain
 
    USE sbc_oce         ! Surface boundary condition: ocean fields
@@ -48,5 +47,4 @@
    USE limvar          ! Ice variables switch
 
-   USE limmsh          ! LIM mesh
    USE limistate       ! LIM initial state
    USE limthd_sal      ! LIM ice thermodynamics: salinity
@@ -310,5 +308,4 @@
       !                                ! Allocate the ice arrays
       ierr =        ice_alloc        ()      ! ice variables
-      ierr = ierr + dom_ice_alloc    ()      ! domain
       ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
       ierr = ierr + thd_ice_alloc    ()      ! thermodynamics
@@ -331,6 +328,4 @@
       !
       CALL lim_thd_sal_init            ! set ice salinity parameters
-      !
-      CALL lim_msh                     ! ice mesh initialization
       !
       IF( ln_limdyn )   CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation
@@ -661,4 +656,6 @@
       diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp ;
       diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp ;
+
+      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
       
    END SUBROUTINE sbc_lim_diag0
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 6746)
@@ -173,12 +173,11 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
                   ! total intermediate advective trends
-                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj)
                   ! update and guess with monotonic sheme
-                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra   * tmask(ji,jj,jk)
-                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
+                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -410,12 +409,11 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
                   ! total intermediate advective trends
-                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj)
                   ! update and guess with monotonic sheme
-                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra
-                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
+                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -438,8 +436,7 @@
          ! --------------------------------------------------
          ! antidiffusive flux on i and j
-
-
-         DO jk = 1, jpkm1
-
+         !
+         DO jk = 1, jpkm1
+            !
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -572,4 +569,5 @@
    END SUBROUTINE tra_adv_tvd_zts
 
+
    SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6746)
@@ -158,4 +158,5 @@
          ELSE                                         ! No restart or restart not found: Euler forward time stepping
             zfact = 1._wp
+            sbc_tsc(:,:,:) = 0._wp
             sbc_tsc_b(:,:,:) = 0._wp
          ENDIF
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6746)
@@ -76,5 +76,6 @@
       REAL(wp) ::   zchl
       REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zqsr100, zqsr_corr
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3
       !!---------------------------------------------------------------------
@@ -83,5 +84,6 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_alloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_alloc( jpi, jpj,      zqsr100, zqsr_corr )
       CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 )
 
@@ -112,8 +114,10 @@
       !                                        !  --------------------------------------
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
-         ! 1% of qsr to compute euphotic layer
+         !                                       ! 1% of qsr to compute euphotic layer
          zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr
          !
-         CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         !
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
          !
          DO jk = 1, nksrp      
@@ -123,5 +127,7 @@
          END DO
          !
-         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         !
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
          !
          DO jk = 1, nksrp      
@@ -133,5 +139,7 @@
          zqsr100(:,:) = 0.01 * qsr(:,:)
          !
-         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         !
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
          !
          DO jk = 1, nksrp      
@@ -226,5 +234,6 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_dealloc( jpi, jpj,      zqsr100, zqsr_corr )
       CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 )
       !
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 6746)
@@ -136,5 +136,5 @@
                   zval = MAX( 1., zstrn(ji,jj) )
                   zval = 1.5 * zval / ( 12. + zval )
-                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
+                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval * ( 1. - fr_i(ji,jj) )
                   zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
                ENDIF
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 6746)
@@ -35,5 +35,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   restotr   ! restoring coeff. on tracers (s-1)
 
-   INTEGER, PARAMETER           ::   npncts   = 5        ! number of closed sea
+   INTEGER, PARAMETER           ::   npncts   = 8        ! number of closed sea
    INTEGER, DIMENSION(npncts)   ::   nctsi1, nctsj1      ! south-west closed sea limits (i,j)
    INTEGER, DIMENSION(npncts)   ::   nctsi2, nctsj2      ! north-east closed sea limits (i,j)
@@ -107,6 +107,5 @@
                
                jl = n_trc_index(jn) 
-               CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000
-               ztrcdta(:,:,:) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl)
+               CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
 
                SELECT CASE ( nn_zdmp_tr )
@@ -187,4 +186,5 @@
       INTEGER :: ji , jj, jk, jn, jl, jc                     ! dummy loop indicesa
       INTEGER :: isrow                                      ! local index
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta       ! 3D  workspace
 
       !!----------------------------------------------------------------------
@@ -207,6 +207,27 @@
             !
                                                         ! Caspian Sea
-            nctsi1(1)   = 332  ; nctsj1(1)   = 243 - isrow
-            nctsi2(1)   = 344  ; nctsj2(1)   = 275 - isrow
+            nctsi1(1)   = 333  ; nctsj1(1)   = 243 - isrow
+            nctsi2(1)   = 342  ; nctsj2(1)   = 274 - isrow
+            !                                           ! Lake Superior
+            nctsi1(2)   = 198  ; nctsj1(2)   = 258 - isrow
+            nctsi2(2)   = 204  ; nctsj2(2)   = 262 - isrow
+            !                                           ! Lake Michigan
+            nctsi1(3)   = 201  ; nctsj1(3)   = 250 - isrow
+            nctsi2(3)   = 203  ; nctsj2(3)   = 256 - isrow
+            !                                           ! Lake Huron
+            nctsi1(4)   = 204  ; nctsj1(4)   = 252 - isrow
+            nctsi2(4)   = 209  ; nctsj2(4)   = 256 - isrow
+            !                                           ! Lake Erie
+            nctsi1(5)   = 206  ; nctsj1(5)   = 249 - isrow
+            nctsi2(5)   = 209  ; nctsj2(5)   = 251 - isrow
+            !                                           ! Lake Ontario
+            nctsi1(6)   = 210  ; nctsj1(6)   = 252 - isrow
+            nctsi2(6)   = 212  ; nctsj2(6)   = 252 - isrow
+            !                                           ! Victoria Lake
+            nctsi1(7)   = 321  ; nctsj1(7)   = 180 - isrow
+            nctsi2(7)   = 322  ; nctsj2(7)   = 189 - isrow
+            !                                           ! Baltic Sea
+            nctsi1(8)   = 297  ; nctsj1(8)   = 270 - isrow
+            nctsi2(8)   = 308  ; nctsj2(8)   = 293 - isrow
             !                                        
             !                                           ! =======================
@@ -277,13 +298,15 @@
          IF(lwp)  WRITE(numout,*)
          !
+         CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )   ! Memory allocation
+         !
          DO jn = 1, jptra
             IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
                 jl = n_trc_index(jn)
-                CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000
+                CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
                 DO jc = 1, npncts
                    DO jk = 1, jpkm1
                       DO jj = nctsj1(jc), nctsj2(jc)
                          DO ji = nctsi1(jc), nctsi2(jc)
-                            trn(ji,jj,jk,jn) = sf_trcdta(jl)%fnow(ji,jj,jk) * tmask(ji,jj,jk) * rf_trfac(jl)
+                            trn(ji,jj,jk,jn) = ztrcdta(ji,jj,jk)
                             trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
                          ENDDO
@@ -293,5 +316,5 @@
              ENDIF
           ENDDO
-          !
+          CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta )
       ENDIF
       !
@@ -313,4 +336,6 @@
       IF( nn_timing == 1 )  CALL timing_start('trc_dmp_init')
       !
+      !Allocate arrays
+      IF( trc_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_dmp_init: unable to allocate arrays' )
 
       IF( lzoom )   nn_zdmp_tr = 0           ! restoring to climatology at closed north or south boundaries
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcdta.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 6746)
@@ -77,5 +77,5 @@
       ALLOCATE( n_trc_index(ntrc), slf_i(ntrc), STAT=ierr0 )
       IF( ierr0 > 0 ) THEN
-         CALL ctl_stop( 'trc_nam: unable to allocate n_trc_index' )   ;   RETURN
+         CALL ctl_stop( 'trc_dta_init: unable to allocate n_trc_index' )   ;   RETURN
       ENDIF
       nb_trcdta      = 0
@@ -91,4 +91,6 @@
       IF(lwp) THEN
          WRITE(numout,*) ' '
+         WRITE(numout,*) 'trc_dta_init : Passive tracers Initial Conditions '
+         WRITE(numout,*) '~~~~~~~~~~~~~~ '
          WRITE(numout,*) ' number of passive tracers to be initialize by data :', ntra
          WRITE(numout,*) ' '
@@ -107,14 +109,18 @@
          DO jn = 1, ntrc
             IF( ln_trc_ini(jn) )  THEN    ! open input file only if ln_trc_ini(jn) is true
-               clndta = TRIM( sn_trcdta(jn)%clvar ) 
-               clntrc = TRIM( ctrcnm   (jn)       ) 
+               clndta = TRIM( sn_trcdta(jn)%clvar )
+               if (jn > jptra) then
+                  clntrc='Dummy' ! By pass weird formats in ocean.output if ntrc > jptra
+               else
+                  clntrc = TRIM( ctrcnm   (jn)       )
+               endif
                zfact  = rn_trfac(jn)
-               IF( clndta /=  clntrc ) THEN 
-                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation :  ',   &
-                  &              'the variable name in the data file : '//clndta//   & 
-                  &              '  must be the same than the name of the passive tracer : '//clntrc//' ')
+               IF( clndta /=  clntrc ) THEN
+                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation    ',   &
+                  &              'Input name of data file : '//TRIM(clndta)//   &
+                  &              ' differs from that of tracer : '//TRIM(clntrc)//' ')
                ENDIF
-               WRITE(numout,*) ' read an initial file for passive tracer number :', jn, ' name : ', clndta, & 
-               &               ' multiplicative factor : ', zfact
+               WRITE(numout,'(a, i4,3a,e11.3)') ' Read IC file for tracer number :', &
+               &            jn, ', name : ', TRIM(clndta), ', Multiplicative Scaling factor : ', zfact
             ENDIF
          END DO
@@ -124,5 +130,5 @@
          ALLOCATE( sf_trcdta(nb_trcdta), rf_trfac(nb_trcdta), STAT=ierr1 )
          IF( ierr1 > 0 ) THEN
-            CALL ctl_stop( 'trc_dta_ini: unable to allocate  sf_trcdta structure' )   ;   RETURN
+            CALL ctl_stop( 'trc_dta_init: unable to allocate  sf_trcdta structure' )   ;   RETURN
          ENDIF
          !
@@ -135,5 +141,5 @@
                IF( sn_trcdta(jn)%ln_tint )  ALLOCATE( sf_trcdta(jl)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
                IF( ierr2 + ierr3 > 0 ) THEN
-                 CALL ctl_stop( 'trc_dta : unable to allocate passive tracer data arrays' )   ;   RETURN
+                 CALL ctl_stop( 'trc_dta_init : unable to allocate passive tracer data arrays' )   ;   RETURN
                ENDIF
             ENDIF
@@ -141,5 +147,5 @@
          ENDDO
          !                         ! fill sf_trcdta with slf_i and control print
-         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta', 'Passive tracer data', 'namtrc' )
+         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta_init', 'Passive tracer data', 'namtrc' )
          !
       ENDIF
@@ -151,5 +157,5 @@
 
 
-   SUBROUTINE trc_dta( kt, sf_dta )
+   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc)
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_dta  ***
@@ -164,9 +170,12 @@
       !!----------------------------------------------------------------------
       INTEGER                     , INTENT(in   ) ::   kt     ! ocean time-step
-      TYPE(FLD), DIMENSION(1)   , INTENT(inout) ::   sf_dta     ! array of information on the field to read
+      TYPE(FLD), DIMENSION(1)     , INTENT(inout) ::   sf_dta     ! array of information on the field to read
+      REAL(wp)                    , INTENT(in   ) ::   ptrfac  ! multiplication factor
+      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL  , INTENT(out  ) ::   ptrc
       !
       INTEGER ::   ji, jj, jk, jl, jkk, ik    ! dummy loop indices
       REAL(wp)::   zl, zi
       REAL(wp), DIMENSION(jpk) ::  ztp                ! 1D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta   ! 3D  workspace
       CHARACTER(len=100) :: clndta
       !!----------------------------------------------------------------------
@@ -176,5 +185,8 @@
       IF( nb_trcdta > 0 ) THEN
          !
+         CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )    ! Memory allocation
+         !
          CALL fld_read( kt, 1, sf_dta )      !==   read data at kt time step   ==!
+         ztrcdta(:,:,:) = sf_dta(1)%fnow(:,:,:) * tmask(:,:,:)    ! Mask
          !
          IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
@@ -185,51 +197,58 @@
             ENDIF
             !
-               DO jj = 1, jpj                         ! vertical interpolation of T & S
+            DO jj = 1, jpj                         ! vertical interpolation of T & S
+               DO ji = 1, jpi
+                  DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+                     zl = fsdept_n(ji,jj,jk)
+                     IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
+                        ztp(jk) = ztrcdta(ji,jj,1)
+                     ELSEIF( zl > gdept_1d(jpk) ) THEN         ! below the last level of data
+                        ztp(jk) =  ztrcdta(ji,jj,jpkm1)
+                     ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
+                        DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
+                           IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
+                              zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
+                              ztp(jk) = ztrcdta(ji,jj,jkk) + ( ztrcdta(ji,jj,jkk+1) - &
+                                        ztrcdta(ji,jj,jkk) ) * zi 
+                           ENDIF
+                        END DO
+                     ENDIF
+                  END DO
+                  DO jk = 1, jpkm1
+                    ztrcdta(ji,jj,jk) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
+                  END DO
+                  ztrcdta(ji,jj,jpk) = 0._wp
+                END DO
+            END DO
+            ! 
+         ELSE                                !==   z- or zps- coordinate   ==!
+            !
+            IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
+               DO jj = 1, jpj
                   DO ji = 1, jpi
-                     DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                        zl = fsdept_n(ji,jj,jk)
-                        IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
-                           ztp(jk) =  sf_dta(1)%fnow(ji,jj,1)
-                        ELSEIF( zl > gdept_1d(jpk) ) THEN         ! below the last level of data
-                           ztp(jk) =  sf_dta(1)%fnow(ji,jj,jpkm1)
-                        ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                           DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                              IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                                 zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                                 ztp(jk) = sf_dta(1)%fnow(ji,jj,jkk) + ( sf_dta(1)%fnow(ji,jj,jkk+1) - &
-                                           sf_dta(1)%fnow(ji,jj,jkk) ) * zi 
-                              ENDIF
-                           END DO
-                        ENDIF
-                     END DO
-                     DO jk = 1, jpkm1
-                        sf_dta(1)%fnow(ji,jj,jk) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                     END DO
-                     sf_dta(1)%fnow(ji,jj,jpk) = 0._wp
+                     ik = mbkt(ji,jj) 
+                     IF( ik > 1 ) THEN
+                        zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik-1)
+                     ENDIF
+                     ik = mikt(ji,jj)
+                     IF( ik > 1 ) THEN
+                        zl = ( fsdept_n(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
+                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik+1)
+                     ENDIF
                   END DO
                END DO
-            ! 
-         ELSE                                !==   z- or zps- coordinate   ==!
-            !                             
-               sf_dta(1)%fnow(:,:,:) = sf_dta(1)%fnow(:,:,:) * tmask(:,:,:)    ! Mask
-               !
-               IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        ik = mbkt(ji,jj) 
-                        IF( ik > 1 ) THEN
-                           zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                           sf_dta(1)%fnow(ji,jj,ik) = (1.-zl) * sf_dta(1)%fnow(ji,jj,ik) + zl * sf_dta(1)%fnow(ji,jj,ik-1)
-                        ENDIF
-                        ik = mikt(ji,jj)
-                        IF( ik > 1 ) THEN
-                           zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
-                           sf_dta(1)%fnow(ji,jj,ik) = (1.-zl) * sf_dta(1)%fnow(ji,jj,ik) + zl * sf_dta(1)%fnow(ji,jj,ik+1)
-                        ENDIF
-                     END DO
-                  END DO
-               ENDIF
-            !
-         ENDIF
+            ENDIF
+            !
+         ENDIF
+         !
+         ! Add multiplicative factor
+         ztrcdta(:,:,:) = ztrcdta(:,:,:) * ptrfac
+         !
+         ! Data structure for trc_ini (and BFMv5.1 coupling)
+         IF( .NOT. PRESENT(ptrc) ) sf_dta(1)%fnow(:,:,:) = ztrcdta(:,:,:)
+         !
+         ! Data structure for trc_dmp
+         IF( PRESENT(ptrc) )  ptrc(:,:,:) = ztrcdta(:,:,:)
          !
          IF( lwp .AND. kt == nit000 ) THEN
@@ -238,11 +257,14 @@
                WRITE(numout,*)
                WRITE(numout,*)'  level = 1'
-               CALL prihre( sf_dta(1)%fnow(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+               CALL prihre( ztrcdta(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
                WRITE(numout,*)'  level = ', jpk/2
-               CALL prihre( sf_dta(1)%fnow(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+               CALL prihre( ztrcdta(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
                WRITE(numout,*)'  level = ', jpkm1
-               CALL prihre( sf_dta(1)%fnow(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+               CALL prihre( ztrcdta(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
                WRITE(numout,*)
          ENDIF
+         !
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta )
+         !
       ENDIF
       !
@@ -255,5 +277,5 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_dta( kt, sf_dta, zrf_trfac )        ! Empty routine
+   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc)        ! Empty routine
       WRITE(*,*) 'trc_dta: You should not have seen this print! error?', kt
    END SUBROUTINE trc_dta
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 6745)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 6746)
@@ -123,7 +123,6 @@
                IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
                   jl = n_trc_index(jn) 
-                  CALL trc_dta( nit000, sf_trcdta(jl) )   ! read tracer data at nit000
-                  trn(:,:,:,jn) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl)
-                  !
+                  CALL trc_dta( nit000, sf_trcdta(jl), rf_trfac(jl) )   ! read tracer data at nit000
+                  trn(:,:,:,jn) = sf_trcdta(jl)%fnow(:,:,:) 
                   IF( .NOT.ln_trcdmp .AND. .NOT.ln_trcdmp_clo ) THEN      !== deallocate data structure   ==!
                      !                                                    (data used only for initialisation)
