Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6963)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6964)
@@ -10,4 +10,5 @@
    !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
    !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
+   !!            3.6  !  2016-06  (C. Rousset) Rewriting (conserves energy)
    !!----------------------------------------------------------------------
 #if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
@@ -95,9 +96,8 @@
       !!                 coriolis terms of the momentum equation
       !!              3) Solve the momentum equation (iterative procedure)
-      !!              4) Prevent high velocities if the ice is thin
-      !!              5) Recompute invariants of the strain rate tensor
+      !!              4) Recompute invariants of the strain rate tensor
       !!                 which are inputs of the ITD, store stress
       !!                 for the next time step
-      !!              6) Control prints of residual (convergence)
+      !!              5) Control prints of residual (convergence)
       !!                 and charge ellipse.
       !!                 The user should make sure that the parameters
@@ -106,4 +106,7 @@
       !!                 e.g. in the Canadian Archipelago
       !!
+      !! ** Notes   : Boundary condition for ice is chosen no-slip 
+      !!              but can be adjusted with param rn_shlat
+      !!
       !! References : Hunke and Dukowicz, JPO97
       !!              Bouillon et al., Ocean Modelling 2009
@@ -115,49 +118,48 @@
       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) ::   zdtevp, 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, zmassU, zmassV                            ! ice/snow mass
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
+      REAL(wp) ::   zTauO, zTauE, zCor                                       ! 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(:,:) ::   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(:,:) ::   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(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses
+      
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! 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(:,:) ::   zswitchU, zswitchV              ! dummy arrays
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice
+
+      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
+      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
+      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
       !!-------------------------------------------------------------------
 
-      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                 )
+      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
+      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
+      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
+      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
 
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -176,59 +178,88 @@
       !
       !------------------------------------------------------------------------------!
-      ! 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
-
+      ! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj) 
+      !------------------------------------------------------------------------------!
+      ! ocean/land mask
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1      ! NO vector opt.
+            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+         END DO
+      END DO
+      CALL lbc_lnk( zfmask, 'F', 1._wp )
+
+      ! Lateral boundary conditions on velocity (modify zfmask)
+      zwf(:,:) = zfmask(:,:)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            IF( zfmask(ji,jj) == 0._wp ) THEN
+               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
+            ENDIF
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         IF( zfmask(1,jj) == 0._wp ) THEN
+            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+         ENDIF
+         IF( zfmask(jpi,jj) == 0._wp ) THEN
+            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
+         ENDIF
+      END DO
+      DO ji = 2, jpim1
+         IF( zfmask(ji,1) == 0._wp ) THEN
+            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+         ENDIF
+         IF( zfmask(ji,jpj) == 0._wp ) THEN
+            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+         ENDIF
+      END DO
+      CALL lbc_lnk( zfmask, 'F', 1._wp )
+
+      !------------------------------------------------------------------------------!
+      ! 1) define some variables and initialize arrays
+      !------------------------------------------------------------------------------!
+      ! ecc2: square of yield ellipse eccenticrity
+      ecc2    = rn_ecc * rn_ecc
+      z1_ecc2 = 1._wp / ecc2
+
+      ! Time step for subcycling
+      zdtevp   = rdt_ice / REAL( nn_nevp )
+      z1_dtevp = 1._wp / zdtevp
+
+      ! alpha parameters (Bouillon 2009)
 #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
+      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+#else
+      zalph1 = ( 2._wp * telast ) * z1_dtevp
+#endif
+      zalph2 = zalph1 * z1_ecc2
+
+      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+
+      ! Initialise stress tensor 
+      zs1 (:,:) = stress1_i (:,:) 
+      zs2 (:,:) = stress2_i (:,:)
+      zs12(:,:) = stress12_i(:,:)
+
+      ! Ice strength
 #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 lim_itd_me_icestrength( nn_icestr )
+      zpresh(:,:) = tmask(:,:,1) *  strength(:,:)
+#else
+      zpresh(:,:) = tmask(:,:,1) *  pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) )
+#endif
+
+      ! scale factors
+      DO jj = k_j1+1, k_jpj-1
+         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 ==!
@@ -242,5 +273,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 ==!
@@ -251,53 +282,48 @@
          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
+            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
+            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(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) )
+            zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
+            zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ocean currents at U-V points
+            v_oceU(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_oceV(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 T points (m*f)
+            zmf(ji,jj)      = zm1 * fcor(ji,jj)
+
+            ! m/dt
+            zmU_t(ji,jj)    = zmassU * z1_dtevp
+            zmV_t(ji,jj)    = zmassV * z1_dtevp
+
+            ! Drag ice-atm.
+            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+
+            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
+            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! masks
+            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
+            ! switches
+            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
+            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
 
          END DO
       END DO
-
+      CALL lbc_lnk( zmf, 'T', 1. )
       !
       !------------------------------------------------------------------------------!
@@ -305,246 +331,240 @@
       !------------------------------------------------------------------------------!
       !
-      ! 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)   &
-                  &            ) * 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)
-
-               !
+         IF(ln_ctl) THEN   ! Convergence test
+            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
+         ENDIF
+
+         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+         DO jj = k_j1, k_jpj-1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+            DO ji = 1, jpim1
+
+               ! 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) * zfmask(ji,jj)
+
+            END DO
+         END DO
+         CALL lbc_lnk( zds, 'F', 1. )
+
          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. )
+            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)
+              
+               ! divergence at T points
+               zdiv  = ( 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)
+               zdiv2 = zdiv * zdiv
+               
+               ! tension at T points
+               zdt  = ( ( 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)
+               zdt2 = zdt * zdt
+               
+               ! delta at T points
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )  
+
+               ! P/delta at T points
+               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
+               
+               ! stress at T points
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
+             
+            END DO
+         END DO
+         CALL lbc_lnk( zp_delt, 'T', 1. )
+
+         DO jj = k_j1, k_jpj-1
+            DO ji = 1, jpim1
+
+               ! 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 at 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)
+         ! --- 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)
+            DO ji = fs_2, fs_jpim1               
+
+               ! U points
+               zfU(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
+               zfV(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)
+
+               ! u_ice at V point
+               u_iceV(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_iceU(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
          !
-         ! Computation of ice velocity
-         !
-         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
-         !
-         IF (MOD(jter,2).eq.0) THEN 
-
+         ! --- Computation of ice velocity --- !
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
+         !  Bouillon et al. 2009 (eq 34-35) => stable
+         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
+            
             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 
+
+                  ! tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+
+                  ! Coriolis at V-points (energy conserving formulation)
+                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
                END DO
             END DO
-
-            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
+            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_bdy
+            CALL bdy_ice_lim_dyn( 'V' )
+#endif         
+
+            DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+                               
+                  ! tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+
+                  ! Coriolis at U-points (energy conserving formulation)
+                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
+                     &           ) * zmaskU(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_bdy
-         CALL bdy_ice_lim_dyn( 'U' )
+            CALL bdy_ice_lim_dyn( 'U' )
 #endif         
+
+         ELSE ! odd iterations
 
             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
+                               
+                  ! tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+
+                  ! Coriolis at U-points (energy conserving formulation)
+                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
+                     &           ) * zmaskU(ji,jj)
                END DO
             END DO
-
-            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
+            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_bdy
+            CALL bdy_ice_lim_dyn( 'U' )
+#endif         
+
+           DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+
+                  ! tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+
+                  ! Coriolis at V-points (energy conserving formulation)
+                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
+                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+               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_bdy
-         CALL bdy_ice_lim_dyn( 'V' )
+            CALL bdy_ice_lim_dyn( 'V' )
 #endif         
 
-         ELSE 
+         ENDIF
+         
+         IF(ln_ctl) THEN   ! Convergence test
             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_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_bdy
-         CALL bdy_ice_lim_dyn( 'U' )
-#endif         
-
-         ENDIF
-         
-         IF(ln_ctl) THEN
-            !---  Convergence test.
-            DO jj = k_j1+1 , k_jpj-1
                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
             END DO
@@ -552,5 +572,5 @@
             IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
          ENDIF
-
+         !
          !                                                ! ==================== !
       END DO                                              !  end loop over jter  !
@@ -558,92 +578,50 @@
       !
       !------------------------------------------------------------------------------!
-      ! 4) Prevent ice velocities when the 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 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
+      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !------------------------------------------------------------------------------!
+      DO jj = k_j1, k_jpj-1 
+         DO ji = 1, jpim1
+
+            ! 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) * zfmask(ji,jj)
+
+         END DO
+      END DO           
+      CALL lbc_lnk( zds, 'F', 1. )
+      
+      DO jj = k_j1+1, k_jpj-1 
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! tension**2 at T points
+            zdt  = ( ( 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)
+            zdt2 = zdt * zdt
+            
+            ! 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)
+            
+            ! shear at T points
+            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! 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)
+            
+            ! delta at T points
+            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )  
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
+
          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_bdy
-      CALL bdy_ice_lim_dyn( 'U' )
-      CALL bdy_ice_lim_dyn( 'V' )
-#endif         
-
-      DO jj = k_j1+1, k_jpj-1 
-         DO ji = fs_2, fs_jpim1
-            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 
-         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
-         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( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
       stress1_i (:,:) = zs1 (:,:)
       stress2_i (:,:) = zs2 (:,:)
@@ -652,5 +630,5 @@
       !
       !------------------------------------------------------------------------------!
-      ! 6) Control prints of residual and charge ellipse
+      ! 5) Control prints of residual and charge ellipse
       !------------------------------------------------------------------------------!
       !
@@ -675,6 +653,6 @@
                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) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
+                     zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 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)
@@ -687,8 +665,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, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
+      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
+      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
+      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
 
    END SUBROUTINE lim_rhg
