Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12740)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12741)
@@ -173,4 +173,7 @@
           <field id="frq_m"    unit="-"    />
 
+          <!-- rheology convergence tests -->
+          <field id="uice_cvg"   long_name="sea ice velocity convergence"   standard_name="sea_ice_velocity_convergence"   unit="m/s" />
+
 	  <!-- ================= -->
           <!-- Add-ons for SIMIP -->
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12740)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12741)
@@ -97,4 +97,5 @@
       rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
                                       !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+   ln_rhg_chkcvg    = .false.         !  check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg)
 /
 !------------------------------------------------------------------------------
@@ -249,5 +250,6 @@
    ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
    ln_icectl        = .false.         !  ice points output for debug (T or F)
-   iiceprt          =  10             !  i-index for debug
-   jiceprt          =  10             !  j-index for debug
-/
+      iiceprt       =  10             !     i-index for debug
+      jiceprt       =  10             !     j-index for debug
+   ln_icediacvg     = .false.         !  check convergence of ice rheology (outputs: file ice_cvg.nc & variable uice_cvg)
+/
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdia
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdia	(revision 12740)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdia	(revision 12741)
@@ -8,5 +8,5 @@
    ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
    ln_icectl        = .false.         !  ice points output for debug (T or F)
-   iiceprt          =  10             !  i-index for debug
-   jiceprt          =  10             !  j-index for debug
+      iiceprt       =  10             !     i-index for debug
+      jiceprt       =  10             !     j-index for debug
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdyn_rhg
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdyn_rhg	(revision 12740)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdyn_rhg	(revision 12741)
@@ -9,3 +9,4 @@
       rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
                                       !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+   ln_rhg_chkcvg    = .false.         !  check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg)
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12740)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12741)
@@ -155,4 +155,5 @@
    INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
+   LOGICAL , PUBLIC ::   ln_rhg_chkcvg    !: check ice rheology convergence 
    !
    !                                     !!** ice-advection namelist (namdyn_adv) **
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg.F90	(revision 12740)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg.F90	(revision 12741)
@@ -110,5 +110,5 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg
       !!-------------------------------------------------------------------
       !
@@ -126,10 +126,11 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          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
-         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp    = ', nn_nevp
-         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast  = ', rn_relast
+         WRITE(numout,*) '      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
+         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp       = ', nn_nevp
+         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast     = ', rn_relast
+         WRITE(numout,*) '      check convergence of rheology                        ln_rhg_chkcvg = ', ln_rhg_chkcvg
       ENDIF
       !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90	(revision 12740)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90	(revision 12741)
@@ -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)
@@ -125,5 +131,4 @@
       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
@@ -141,5 +146,4 @@
       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
@@ -160,6 +164,7 @@
       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
       !! --- SIMIP diags
@@ -174,5 +179,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
@@ -345,10 +359,13 @@
          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(ln_rhg_chkcvg) 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) --- !
@@ -667,11 +684,8 @@
          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(ln_rhg_chkcvg) THEN
+            CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
+         ENDIF
          !
          !                                                ! ==================== !
@@ -734,10 +748,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. &
@@ -796,5 +804,5 @@
          DEALLOCATE( zsig1 , zsig2 , zsig3 )
       ENDIF
-      
+
       ! --- SIMIP --- !
       IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
@@ -852,5 +860,85 @@
       ENDIF
       !
+      ! --- convergence tests --- !
+      IF( ln_rhg_chkcvg ) THEN
+         IF( iom_use('uice_cvg') ) THEN  ! output: u(t=nn_nevp) - u(t=nn_nevp-1)
+            CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
+               &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
+         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 )   istatus = NF90_CLOSE(ncvgid)
+      ENDIF
+      
+   END SUBROUTINE rhg_cvg
 
 
@@ -910,4 +998,5 @@
    END SUBROUTINE rhg_evp_rst
 
+   
 #else
    !!----------------------------------------------------------------------
