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 6789)
+++ branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6796)
@@ -10,5 +10,5 @@
    !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
    !!            3.5  !  2012-08  (R. Benshila)  AGRIF
-   !!            3.6  !  2016-06  (C. Rousset) Rewriting and add the possibility to use mEVP from Bouillon 2013
+   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -112,10 +112,10 @@
       CHARACTER (len=50) ::   charout
 
-      REAL(wp) ::   dtevp, z1_dtevp                                          ! time step for subcycling
+      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                                            ! ice/snow mass
+      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, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel   ! temporary scalars
+      REAL(wp) ::   zTauO, zTauB, zTauE, zCor, zvel                          ! temporary scalars
 
       REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
@@ -123,14 +123,14 @@
       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(:,:) ::   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
@@ -148,7 +148,7 @@
       !!-------------------------------------------------------------------
 
-      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, 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, zswitch1, zswitch2, zfmask, zwf )
@@ -205,6 +205,6 @@
 
       ! Time step for subcycling
-      dtevp    = rdt_ice / REAL( nn_nevp )
-      z1_dtevp = 1._wp / dtevp
+      zdtevp   = rdt_ice / REAL( nn_nevp )
+      z1_dtevp = 1._wp / zdtevp
 
       ! alpha parameters (Bouillon 2009)
@@ -228,5 +228,4 @@
       ! Ice strength
       CALL lim_itd_me_icestrength( nn_icestr )
-      zpresh(:,:) = tmask(:,:,1) *  strength(:,:)
 
       ! scale factors
@@ -263,6 +262,6 @@
 
             ! 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)
+            zaU(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)
+            zaV(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
@@ -270,28 +269,36 @@
             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)
+            zmassU = ( zm1 * e1t(ji+1,jj) + zm2 * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+            zmassV = ( 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  ) )
+            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 * ff_t(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) = - 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)
+            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)
 
             ! 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) ) )
+            zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )
+            zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )
 
          END DO
       END DO
+      CALL lbc_lnk( zmf, 'T', 1. )
       !
       !------------------------------------------------------------------------------!
@@ -311,5 +318,5 @@
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
          DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
-            DO ji = 1, fs_jpim1
+            DO ji = 1, jpim1
 
                ! shear at F points
@@ -346,5 +353,5 @@
 
                ! P/delta at T points
-               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
                
                ! stress at T points
@@ -375,5 +382,5 @@
 
                ! U points
-               zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+               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)                                                                      &
@@ -383,5 +390,5 @@
 
                ! V points
-               zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+               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)                                                                      &
@@ -391,9 +398,9 @@
 
                ! 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_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_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+               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)
 
@@ -402,5 +409,5 @@
          !
          ! --- 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. 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
@@ -409,20 +416,33 @@
                DO ji = fs_2, fs_jpim1
 
-                  zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) )
+                  ! 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) ) )
+
+                  ! tau_bottom/v_ice
+                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
+                  ztauB = - tau_icebfr(ji,jj) / zvel
+
+                  ! 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) )
+
+                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * v_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   
-                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj)
-                  zTauO   =   za2(ji,jj) * rhoco * 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) ) )
-                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
-                  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)
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
+                     &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &           ) * zswitch2(ji,jj)
                   ! 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)
+                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
+                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
 
                END DO
@@ -439,25 +459,34 @@
             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) 
+                               
+                  ! 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) ) )
+
+                  ! tau_bottom/u_ice
+                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
+                  ztauB = - tau_icebfr(ji,jj) / zvel
+
+                  ! 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) ) )
                   
-                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
-                  
-                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj)
-                  zTauO   =   za1(ji,jj) * rhoco * 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) ) )
-                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
-                  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)
+                  ! 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) )
+
+                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * u_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
+                     &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &           ) * zswitch1(ji,jj)
                   ! 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)
+                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
+                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
                END DO
             END DO
@@ -476,20 +505,33 @@
                DO ji = fs_2, fs_jpim1
 
-                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
-
-                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj)
-                  zTauO   =   za1(ji,jj) * rhoco * 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) ) )
-                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
-                  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)
+                  ! 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) ) )
+
+                  ! tau_bottom/u_ice
+                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
+                  ztauB = - tau_icebfr(ji,jj) / zvel
+
+                  ! 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) )
+
+                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * u_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
+                     &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast
+                     &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &           ) * zswitch1(ji,jj)
                   ! 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)
+                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
+                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj)
                END DO
             END DO
@@ -505,25 +547,34 @@
             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_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) )
          
-                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj)
-                  zTauO   =   za2(ji,jj) * rhoco * 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) ) )
-                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi )
-                  zCor    =   zcor2(ji,jj)  * zu_ice2
-                  zmt     =   zmass2(ji,jj) * z1_dtevp
+                  ! 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) ) )
+
+                  ! tau_bottom/v_ice
+                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
+                  ztauB = - tau_icebfr(ji,jj) / zvel
                   
-                  ! 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)
+                  ! 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) )
+
+                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * v_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
+                     &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast
+                     &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &           ) * zswitch2(ji,jj)
                   ! 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)
+                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
+                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj)
                END DO
             END DO
@@ -579,5 +630,5 @@
       !------------------------------------------------------------------------------!
       DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1
+         DO ji = 1, jpim1
 
             ! shear at F points
@@ -649,7 +700,7 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  IF (zpresh(ji,jj) > 1.0) THEN
-                     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) )
+                  IF (strength(ji,jj) > 1.0) THEN
+                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
+                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )
                      WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
                      CALL prt_ctl_info(charout)
@@ -662,7 +713,7 @@
       ENDIF
       !
-      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, z1_e1t0, z1_e2t0, zp_delt )
+      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
+      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
       CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
       CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf )
