Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_evp.F90	(revision 11843)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_evp.F90	(revision 13662)
@@ -41,4 +41,5 @@
    USE prtctl         ! Print control
 
+   USE netcdf         ! NetCDF library for convergence test
    IMPLICIT NONE
    PRIVATE
@@ -49,4 +50,9 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+
+   !! for convergence tests
+   INTEGER ::   ncvgid   ! netcdf file id
+   INTEGER ::   nvarid   ! netcdf variable id
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -119,16 +125,17 @@
       REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
       REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
+      REAl(wp) ::   zbetau, zbetav
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
-      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
+      REAL(wp) ::   zp_delf, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars
       REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
       REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
       REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
       !
-      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
       REAL(wp) ::   zintb, zintn                                        ! dummy argument
       REAL(wp) ::   zfac_x, zfac_y
       REAL(wp) ::   zshear, zdum1, zdum2
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
       !
@@ -137,9 +144,8 @@
       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) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
-!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
       !                                                                 !    ocean surface (ssh_m) if ice is not embedded
@@ -155,12 +161,14 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! 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 becomes very small
       REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
+      !! --- check convergence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
       !! --- diags
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
+      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p         
       !! --- SIMIP diags
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
@@ -174,5 +182,14 @@
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
       !
-!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
+      ! for diagnostics and convergence tests
+      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
+            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
+         END DO
+      END DO
+      !
+      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
       !------------------------------------------------------------------------------!
       ! 0) mask at F points for the ice
@@ -187,9 +204,9 @@
 
       ! 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) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
+               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
+                  &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
             ENDIF
          END DO
@@ -197,16 +214,16 @@
       DO jj = 2, jpjm1
          IF( zfmask(1,jj) == 0._wp ) THEN
-            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
          ENDIF
          IF( zfmask(jpi,jj) == 0._wp ) THEN
-            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
-         ENDIF
+            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) )
+        ENDIF
       END DO
       DO ji = 2, jpim1
          IF( zfmask(ji,1) == 0._wp ) THEN
-            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
          ENDIF
          IF( zfmask(ji,jpj) == 0._wp ) THEN
-            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) )
          ENDIF
       END DO
@@ -222,16 +239,17 @@
       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( .NOT. ln_aEVP ) THEN
-         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+         zdtevp   = rdt_ice / REAL( nn_nevp )
+         zalph1 =   2._wp * rn_relast * REAL( nn_nevp )
          zalph2 = zalph1 * z1_ecc2
 
          z1_alph1 = 1._wp / ( zalph1 + 1._wp )
          z1_alph2 = 1._wp / ( zalph2 + 1._wp )
-      ENDIF
+      ELSE
+         zdtevp   = rdt_ice
+         ! zalpha parameters set later on adaptatively
+      ENDIF
+      z1_dtevp = 1._wp / zdtevp
          
       ! Initialise stress tensor 
@@ -244,5 +262,5 @@
 
       ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
-      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
+      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
       ELSE                         ;   zkt = 0._wp
       ENDIF
@@ -315,12 +333,12 @@
                zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
                ! ice-bottom stress at U points
-               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
-               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj)
+               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
                ! ice-bottom stress at V points
-               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
-               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj)
+               ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
                ! ice_bottom stress at T points
-               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
-               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+               zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
             END DO
          END DO
@@ -345,13 +363,16 @@
          l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
          !
-!!$         IF(ln_ctl) THEN   ! Convergence test
-!!$            DO jj = 1, jpjm1
-!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
-!!$               zv_ice(:,jj) = v_ice(:,jj)
-!!$            END DO
-!!$         ENDIF
+         ! convergence test
+         IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
+                  zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
+               END DO
+            END DO
+         ENDIF
 
          ! --- 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 jj = 1, jpjm1
             DO ji = 1, jpim1
 
@@ -363,8 +384,7 @@
             END DO
          END DO
-         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
-
-         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
-            DO ji = 2, jpi ! no vector loop
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 ! no vector loop
 
                ! shear**2 at T points (doc eq. A16)
@@ -386,10 +406,31 @@
                
                ! delta at T points
-               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
-
-               ! P/delta at T points
-               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
-
-               ! alpha & beta for aEVP
+               zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1._wp )
+         
+         ! P/delta at T points
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+            END DO
+         END DO
+
+         DO jj = 2, jpj    ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
+            DO ji = 2, jpi ! no vector loop
+
+               ! divergence at T points (duplication to avoid communications)
+               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_e1e2t(ji,jj)
+               
+               ! tension at T points (duplication to avoid communications)
+               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_e1e2t(ji,jj)
+
+               ! alpha for aEVP
                !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
                !   alpha = beta = sqrt(4*gamma)
@@ -399,22 +440,37 @@
                   zalph2   = zalph1
                   z1_alph2 = z1_alph1
+                  ! explicit:
+                  ! z1_alph1 = 1._wp / zalph1
+                  ! z1_alph2 = 1._wp / zalph1
+                  ! zalph1 = zalph1 - 1._wp
+                  ! zalph2 = zalph1
                ENDIF
                
                ! stress at T points (zkt/=0 if landfast)
-               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
-               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+               zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
              
             END DO
          END DO
-         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
-
+
+         ! Save beta at T-points for further computations
+         IF( ln_aEVP ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+               END DO
+            END DO
+         ENDIF
+         
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
 
-               ! alpha & beta for aEVP
+               ! alpha 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) ) )
+                  zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
                   z1_alph2 = 1._wp / ( zalph2 + 1._wp )
-                  zbeta(ji,jj) = zalph2
+                  ! explicit:
+                  ! z1_alph2 = 1._wp / zalph2
+                  ! zalph2 = zalph2 - 1._wp
                ENDIF
                
@@ -486,17 +542,20 @@
                   !
                   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
-                        &                                  + zRHS + 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
-                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
+                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
+                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetav + 1._wp )                                              &
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                         &           )   * zmsk00y(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
-                        &                                     + zRHS + 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
-                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
-                        &            )   * zmsk00y(ji,jj)
+                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                    + zRHS + 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_lf_relax )         & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )  * zmsk00y(ji,jj)
                   ENDIF
                END DO
@@ -537,17 +596,20 @@
                   !
                   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
-                        &                                  + zRHS + 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
-                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
+                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetau + 1._wp )                                              &
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                         &           )   * zmsk00x(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
-                        &                                     + zRHS + 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
-                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
-                        &            )   * zmsk00x(ji,jj)
+                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                    + zRHS + 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_lf_relax )         & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
                   ENDIF
                END DO
@@ -590,17 +652,20 @@
                   !
                   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
-                        &                                  + zRHS + 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
-                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
+                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetau + 1._wp )                                              &
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                         &           )   * zmsk00x(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
-                        &                                     + zRHS + 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
-                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
-                        &            )   * zmsk00x(ji,jj)
+                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                    + zRHS + 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_lf_relax )         & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00x(ji,jj)
                   ENDIF
                END DO
@@ -641,17 +706,20 @@
                   !
                   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
-                        &                                  + zRHS + 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
-                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
+                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
+                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetav + 1._wp )                                              & 
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                         &           )   * zmsk00y(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
-                        &                                     + zRHS + 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
-                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
-                        &            )   * zmsk00y(ji,jj)
+                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                    + zRHS + 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_lf_relax )         & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00y(ji,jj)
                   ENDIF
                END DO
@@ -667,15 +735,11 @@
          ENDIF
 
-!!$         IF(ln_ctl) THEN   ! Convergence test
-!!$            DO jj = 2 , jpjm1
-!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
-!!$            END DO
-!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
-!!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
-!!$         ENDIF
+         ! convergence test
+         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
          !
          !                                                ! ==================== !
       END DO                                              !  end loop over jter  !
       !                                                   ! ==================== !
+      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
       !
       !------------------------------------------------------------------------------!
@@ -702,4 +766,6 @@
             zdt2 = zdt * zdt
             
+            zten_i(ji,jj) = zdt
+
             ! shear**2 at T points (doc eq. A16)
             zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
@@ -716,14 +782,14 @@
             
             ! delta at T points
-            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
-            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
-            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
+            zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta  
+            rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
+            pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
 
          END DO
       END DO
-      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., &
+         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. )
       
       ! --- Store the stress tensor for the next time step --- !
-      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
       pstress1_i (:,:) = zs1 (:,:)
       pstress2_i (:,:) = zs2 (:,:)
@@ -734,10 +800,4 @@
       ! 5) diagnostics
       !------------------------------------------------------------------------------!
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
-         END DO
-      END DO
-
       ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
       IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
@@ -758,43 +818,76 @@
       IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
 
-      ! --- stress tensor --- !
-      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
-         !
-         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
+      ! --- Stress tensor invariants (SIMIP diags) --- !
+      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
+         !
+         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !         
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
-                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
-                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
-
-               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
-
-               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
-
-!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
-!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
-!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
-!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
-               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
-               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
-               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
-            END DO
-         END DO
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
-         !
-         CALL iom_put( 'isig1' , zsig1 )
-         CALL iom_put( 'isig2' , zsig2 )
-         CALL iom_put( 'isig3' , zsig3 )
-         !
-         ! Stress tensor invariants (normal and shear stress N/m)
-         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
-         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
-
-         DEALLOCATE( zsig1 , zsig2 , zsig3 )
-      ENDIF
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+            
+               ! Ice stresses
+               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
+               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
+               ! I know, this can be confusing...
+               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 
+               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
+               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
+               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+               
+               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
+               zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
+               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
+               
+            END DO
+         END DO         
+         !
+         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
+         
+         DEALLOCATE ( zsig_I, zsig_II )
+         
+      ENDIF
+
+      ! --- Normalized stress tensor principal components --- !
+      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
+      ! Recommendation 1 : we use ice strength, not replacement pressure
+      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
+      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
+         !
+         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )         
+         !         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+            
+               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
+               !                        and **deformations** at current iterates
+               !                        following Lemieux & Dupont (2020)
+               zfac             =   zp_delt(ji,jj)
+               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
+               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
+               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+               
+               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
+               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
+               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
       
+               ! Normalized  principal stresses (used to display the ellipse)
+               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
+               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
+            END DO
+         END DO               
+         !
+         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
+         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
+      
+         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
+         
+      ENDIF
+
       ! --- SIMIP --- !
       IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
@@ -852,5 +945,90 @@
       ENDIF
       !
+      ! --- convergence tests --- !
+      IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
+         IF( iom_use('uice_cvg') ) THEN
+            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
+               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
+                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
+            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
+               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
+                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
+            ENDIF
+         ENDIF
+      ENDIF      
+      !
+      DEALLOCATE( zmsk00, zmsk15 )
+      !
    END SUBROUTINE ice_dyn_rhg_evp
+
+
+   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE rhg_cvg  ***
+      !!                     
+      !! ** Purpose :   check convergence of oce rheology
+      !!
+      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
+      !!                during the sub timestepping of rheology so as:
+      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
+      !!                This routine is called every sub-iteration, so it is cpu expensive
+      !!
+      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
+      !!----------------------------------------------------------------------
+      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
+      !!
+      INTEGER           ::   it, idtime, istatus
+      INTEGER           ::   ji, jj          ! dummy loop indices
+      REAL(wp)          ::   zresm           ! local real 
+      CHARACTER(len=20) ::   clname
+      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
+      !!----------------------------------------------------------------------
+
+      ! create file
+      IF( kt == nit000 .AND. kiter == 1 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !
+         IF( lwm ) THEN
+            clname = 'ice_cvg.nc'
+            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
+            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
+            istatus = NF90_ENDDEF(ncvgid)
+         ENDIF
+         !
+      ENDIF
+
+      ! time
+      it = ( kt - 1 ) * kitermax + kiter
+      
+      ! convergence
+      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
+         zresm = 0._wp
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
+                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)
+            END DO
+         END DO
+         zresm = MAXVAL( zres )
+         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
+      ENDIF
+
+      IF( lwm ) THEN
+         ! write variables
+         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
+         ! close file
+         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid)
+      ENDIF
+      
+   END SUBROUTINE rhg_cvg
 
 
@@ -910,4 +1088,5 @@
    END SUBROUTINE rhg_evp_rst
 
+   
 #else
    !!----------------------------------------------------------------------
