Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8626)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8678)
@@ -90,5 +90,6 @@
 !------------------------------------------------------------------------------
    ln_rhg_EVP       = .true.          !  EVP rheology
-      rn_creepl     =   1.0e-12       !     creep limit [1/s]
+      ln_aEVP       = .false.         !     adaptive rheology (Kimmritz et al. 2016 & 2017)
+      rn_creepl     =   2.0e-9        !     creep limit [1/s]
       rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve          
       nn_nevp       = 120             !     number of EVP subcycles                             
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8626)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8678)
@@ -178,4 +178,5 @@
    !
    !                                     !!** ice-rheology namelist (namrhg) **
+   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
    REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90	(revision 8626)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90	(revision 8678)
@@ -78,5 +78,5 @@
       CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic !
          !                             !------------------------!
-         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
          !         
       END SELECT
@@ -107,5 +107,5 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
       !!-------------------------------------------------------------------
       !
@@ -125,4 +125,5 @@
          WRITE(numout,*) '   Namelist namdyn_rhg:'
          WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP = ', ln_rhg_EVP
+         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP    = ', ln_aEVP
          WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl
          WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90	(revision 8626)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90	(revision 8678)
@@ -53,5 +53,5 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, u_ice, v_ice, pshear_i, pdivu_i, pdelta_i )
+   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
@@ -95,15 +95,16 @@
       !!              5) Diagnostics including charge ellipse
       !!
-      !! ** 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)
+      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
+      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
+      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
+      !!              (i.e. more stable and better convergence)
       !!
       !! References : Hunke and Dukowicz, JPO97
       !!              Bouillon et al., Ocean Modelling 2009
       !!              Bouillon et al., Ocean Modelling 2013
+      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! time step
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   u_ice, v_ice 
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i, pdivu_i, pdelta_i 
       !!
@@ -114,5 +115,5 @@
       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) ::   zalph1, z1_alph1, zalph2, z1_alph2                       ! alpha coef from Bouillon 2009 or Kimmritz 2017
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
@@ -126,7 +127,9 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
-      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
@@ -231,16 +234,12 @@
 
       ! 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 )
-
+      IF( .NOT. ln_aEVP ) THEN
+         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+         zalph2 = zalph1 * z1_ecc2
+
+         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+      ENDIF
+         
       ! Initialise stress tensor 
       zs1 (:,:) = pstress1_i (:,:) 
@@ -304,8 +303,11 @@
             zmf(ji,jj)      = zm1 * ff_t(ji,jj)
 
+            ! dt/m at T points (for alpha and beta coefficients)
+            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
+            
             ! 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)
@@ -326,5 +328,5 @@
          END DO
       END DO
-      CALL lbc_lnk( zmf, 'T', 1. )
+      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
       !
       !------------------------------------------------------------------------------!
@@ -380,4 +382,14 @@
                ! P/delta at T points
                zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+
+               ! alpha & beta for aEVP
+               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
+               !   alpha = beta = sqrt(4*gamma)
+               IF( ln_aEVP ) THEN
+                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+                  zalph2   = zalph1
+                  z1_alph2 = z1_alph1
+               ENDIF
                
                ! stress at T points
@@ -392,4 +404,11 @@
             DO ji = 1, jpim1
 
+               ! alpha & beta for aEVP
+               IF( ln_aEVP ) THEN
+                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+                  zbeta(ji,jj) = zalph2
+               ENDIF
+               
                ! 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) )
@@ -411,6 +430,6 @@
                   &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
                   &                  ) * r1_e1e2u(ji,jj)
-                  !
-                  !                !--- V points
+               !
+               !                !--- 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)    &
@@ -419,18 +438,17 @@
                   &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
                   &                  ) * r1_e1e2v(ji,jj)
-                  !
-                  !                !--- u_ice at V point
+               !
+               !                !--- 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_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 --- !
-         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
          !  Bouillon et al. 2009 (eq 34-35) => stable
          IF( MOD(jter,2) == 0 ) THEN ! even iterations
@@ -438,26 +456,33 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                     !                 !--- tau_io/(v_oce - v_ice)
+                  !                 !--- tau_io/(v_oce - v_ice)
                   zTauO = zaV(ji,jj) * zrhoco * 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) ) )
-                     !                 !--- Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-                     !
-                     !                 !--- tau_bottom/v_ice
+                  !
+                  !                 !--- 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)
+                  !
+                  !                 !--- Coriolis at V-points (energy conserving formulation)
                   zCory(ji,jj)  = - 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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-                     !
-                     !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                     !
-                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(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) * ( zbeta(ji,jj) + 1._wp ) + 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
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ELSE               !--- ice velocity using EVP 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)
@@ -466,10 +491,5 @@
                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                      &            ) * zmaskV(ji,jj)
-                     !
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
-                  !
+                  ENDIF
                END DO
             END DO
@@ -483,29 +503,34 @@
             !
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                               
-                  ! tau_io/(u_oce - u_ice)
+               DO ji = fs_2, fs_jpim1          
+                  !                 !--- tau_io/(u_oce - u_ice)
                   zTauO = zaU(ji,jj) * zrhoco * 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) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
+                  !
+                  !                 !--- 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)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(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) * ( zbeta(ji,jj) + 1._wp ) + 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
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &            ) * zmaskU(ji,jj)
+                  ELSE               !--- ice velocity using EVP 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)
@@ -514,8 +539,5 @@
                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                      &            ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -532,38 +554,40 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-
-                  ! tau_io/(u_oce - u_ice)
+                  !                 !--- tau_io/(u_oce - u_ice)
                   zTauO = zaU(ji,jj) * zrhoco * 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) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
+                  !
+                  !                 !--- 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)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(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) * ( zbeta(ji,jj) + 1._wp ) + 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
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &            ) * zmaskU(ji,jj)
+                  ELSE               !--- ice velocity using EVP 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
-                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                      &            ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -578,28 +602,33 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-         
-                  ! tau_io/(v_oce - v_ice)
+                  !                 !--- tau_io/(v_oce - v_ice)
                   zTauO = zaV(ji,jj) * zrhoco * 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) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-
-                  ! tau_bottom/v_ice
+                  !
+                  !                 !--- 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)
+                  !
+                  !                 !--- Coriolis at v-points (energy conserving formulation)
                   zCory(ji,jj)  = - 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
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                  
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(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) * ( zbeta(ji,jj) + 1._wp ) + 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
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ELSE               !--- ice velocity using EVP 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)
@@ -608,8 +637,5 @@
                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                      &            ) * zmaskV(ji,jj)
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
+                  ENDIF
                END DO
             END DO
