Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90	(revision 13522)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90	(revision 13536)
@@ -160,5 +160,4 @@
    INTEGER , PUBLIC ::   nn_nout_vp       !: Number of outer iterations
    INTEGER , PUBLIC ::   nn_ninn_vp       !: Number of inner iterations (linear system solver)
-   LOGICAL , PUBLIC ::   ln_smooth_vp     !: zeta viscosity smoothing (if yes, tanh function of Lemieux and Tremblay JGR 2009)
    LOGICAL , PUBLIC ::   ln_zebra_vp      !: activate zebra (solve the linear system problem every odd j-band, then one every even one)
    REAL(wp), PUBLIC ::   rn_relaxu_vp     !: U-relaxation factor (MV: can probably be merged with V-factor once ok)
Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90	(revision 13522)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90	(revision 13536)
@@ -56,5 +56,5 @@
       !!
       !! ** Action  : comupte - ice velocity (u_ice, v_ice)
-      !!                      - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i)
+      !!                      - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i) - EVP only
       !!                      - shear, divergence and delta (shear_i, divu_i, delta_i)
       !!--------------------------------------------------------------------
@@ -88,5 +88,5 @@
       CASE( np_rhgVP  )                ! Viscous-Plastic rheology                    !
          !                             !---------------------------------------------!
-         CALL ice_dyn_rhg_vp ( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
+         CALL ice_dyn_rhg_vp ( kt, shear_i, divu_i, delta_i )
          !         
       END SELECT
@@ -94,5 +94,5 @@
       IF( lrst_ice ) THEN                       !* write EVP fields in the restart file
          IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt )
-         IF( ln_rhg_VP  )   CALL rhg_lsr_rst( 'WRITE', kt )
+         ! MV note: no restart needed for VP as there is no time equation for stress tensor
       ENDIF
       !
@@ -121,5 +121,5 @@
       !!
       NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg, &           !- evp
-    &                       ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_smooth_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp  !-vp
+    &                       ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp  !-vp
       !!-------------------------------------------------------------------
       
@@ -149,5 +149,4 @@
          WRITE(numout,*) '         number of outer iterations                        nn_nout_vp    = ', nn_nout_vp
          WRITE(numout,*) '         number of inner iterations                        nn_ninn_vp    = ', nn_ninn_vp
-         WRITE(numout,*) '         tanh viscosity smooothing                         ln_smooth_vp  = ', ln_smooth_vp
          WRITE(numout,*) '         activate zebra solver                             ln_zebra_vp   = ', ln_zebra_vp
          WRITE(numout,*) '         relaxation factor for u                           rn_relaxu_vp  = ', rn_relaxu_vp
Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90	(revision 13522)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90	(revision 13536)
@@ -137,4 +137,5 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
       !
@@ -168,5 +169,6 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
       !! --- diags
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
+      REAL(wp) ::   zsig1, zsig2, zsig12
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i         
       !! --- SIMIP diags
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
@@ -406,4 +408,6 @@
                ! delta at T points
                zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+               
+               zdeltastar_t(ji,jj) = zdelta + rn_creepl                     ! store delta* at previous iterate (for ellipse computation)
 
                ! P/delta at T points
@@ -695,4 +699,12 @@
                         &           )   * zmsk00y(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                                     !---> MV : comment above is misleading as EVP is a purely explicit method
+                                     !---> MV : I found code below quite unreadable :-)
+                                     !---> MV : two maskings + landfast / no landfast options, I think this is a bit too much
+                                     !---> MV : I would suggest to split at least the landfast / masking steps if possible
+                                     !---> MV : Should also comment why 1% of ocean velocity is assumed - not sure it makes sense btw
+                                     !---> MV : I would say 100% of ocean current makes much more sense for free-drift
+                                     !---> MV : of very low concentration sea ice..
+                                     !---> MV : but maybe it is a numerical trick... in brief it's n
                      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)
@@ -715,5 +727,5 @@
 
          ! convergence test
-         IF( ln_rhg_chkcvg )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
+         IF( ln_rhg_chkcvg )   CALL rhg_cvg_evp( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
          !
          !                                                ! ==================== !
@@ -745,4 +757,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  )  &
@@ -797,41 +811,78 @@
       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 = 2, jpj - 1
+            DO ji = 2, jpi - 1
+            
+               ! 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 ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress
+               
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
+         
+         !
+         ! 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) )         
+         !         
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
+            
+               ! 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) - zdeltastar_t(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), 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 ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress
+
+               ! Normalized  principal stresses (used to display the ellipse)
+               z1_strength      =   1._wp / strength(ji,jj)
+               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  - zsig_II(ji,jj) ) * z1_strength
+               zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+            END DO
+         END DO
+         
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
+               
+         !
+         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
+         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
+
+         DEALLOCATE( zsig1_p , zsig2_p )
+         
+      ENDIF
       ! --- SIMIP --- !
       IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
@@ -907,9 +958,9 @@
 
 
-   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )
+   SUBROUTINE rhg_cvg_evp( kt, kiter, kitermax, pu, pv, pub, pvb )
       !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE rhg_cvg  ***
+      !!                    ***  ROUTINE rhg_cvg_evp  ***
       !!                     
-      !! ** Purpose :   check convergence of oce rheology
+      !! ** Purpose :   check convergence of evp ice rheology
       !!
       !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
@@ -935,6 +986,6 @@
          IF( lwp ) THEN
             WRITE(numout,*)
-            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
-            WRITE(numout,*) '~~~~~~~'
+            WRITE(numout,*) 'rhg_cvg_evp : ice rheology convergence control'
+            WRITE(numout,*) '~~~~~~~~~~~'
          ENDIF
          !
@@ -974,5 +1025,5 @@
       ENDIF
       
-   END SUBROUTINE rhg_cvg
+   END SUBROUTINE rhg_cvg_evp
 
 
Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13522)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13536)
@@ -1,29 +1,22 @@
-! TODOLIST
-!
-! Define all symbols
-! - Do viscosity smoothing with sum (differentiable rheology)
-! - Re-calculate internal force diagnostic (end of the routine)
-! - Do proper masking of output fileds with ice mass (zmsk00 see EVP routine)
-!
-! quality control - compare code to mitGCM
-! quality control - comparing EVP and VP codes
-! 
+! TODOLOIST
 !
 ! Commit
 ! Compile
 !
+! - Re-calculate internal force diagnostic (currently incorrect, see end of the routine)
+! - QC:compare code to mitGCM
+! - QC: comparing EVP and VP codes 
+!
 ! Clarify
 ! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ?
-! --- Swap of independent term in the U-V solvers ?
-! --- Is stress tensor for restart needed ?
-! --- Is stress tensor calculated at the end of the time step
+! --- Swap of independent term in the U-V solvers ? -> JF Lemieux says maybe
+! --- Is stress tensor for restart needed ? No!!!
 !
 ! Test
-! --- Can we add the 15% mask in the convergence criteria ?
 ! --- Try ADI for u-v solver
 !
 ! Add
 ! - Charge ellipse as an output (good one from Lemieux and Dupont)
-! - Bottom drag
+! - Bottom drag (landfast ice)
 ! - Tensile strength
 ! - agrif
@@ -52,5 +45,4 @@
    !!----------------------------------------------------------------------
    !!   ice_dyn_rhg_vp : computes ice velocities from VP rheolog with LSR solvery
-   !!   rhg_vp_rst     : read/write EVP fields in ice restart
    !!----------------------------------------------------------------------
    USE phycst         ! Physical constants
@@ -79,5 +71,4 @@
 
    PUBLIC   ice_dyn_rhg_vp   ! called by icedyn_rhg.F90
-   PUBLIC   rhg_vp_rst       ! called by icedyn_rhg.F90
 
    !! for convergence tests
@@ -99,5 +90,5 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_rhg_vp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
+   SUBROUTINE ice_dyn_rhg_vp( kt, pshear_i, pdivu_i, pdelta_i )
       !!-------------------------------------------------------------------
       !!
@@ -164,5 +155,4 @@
       !!
       INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
       !!
@@ -180,5 +170,5 @@
       REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
       REAL(wp) ::   zm2, zm3, zmassU, zmassV                            ! ice/snow mass and volume
-      REAL(wp) ::   zdeltat, zdeltat_star, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
+      REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars
       REAL(wp) ::   zp_deltastar_f
       REAL(wp) ::   zfac, zfac1, zfac2, zfac3
@@ -192,5 +182,6 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta* at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zef                             ! Viscosity pre-factor at F point
@@ -205,5 +196,4 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
-      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
       !                                                                 !    ocean surface (ssh_m) if ice is not embedded
@@ -225,6 +215,6 @@
 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
      !
-      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) ::   zmsk01x, zmsk01y                ! mask for lots of ice (1), little ice (0)
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence (1), no ice (0)
       !
       REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
@@ -232,12 +222,13 @@
       REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
       !! --- diags
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3         
-      !! --- SIMIP diags
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)      
+      REAL(wp) ::   zsig1, zsig2, zsig12
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f, zten_i   ! stress tensor components
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s, SIMIP)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s, SIMIP)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s, SIMIP)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s, SIMIP)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s, SIMIP)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s, SIMIP)      
                         
       !!----------------------------------------------------------------------------------------------------------------------
@@ -269,10 +260,5 @@
       ecc2    = rn_ecc * rn_ecc
       z1_ecc2 = 1._wp / ecc2
-         
-      ! Initialise stress tensor from previous time step
-      zs1 (:,:) = pstress1_i (:,:) 
-      zs2 (:,:) = pstress2_i (:,:)
-      zs12(:,:) = pstress12_i(:,:)
-      
+               
       ! Initialise convergence checks
       IF( ln_rhg_chkcvg ) THEN
@@ -389,13 +375,13 @@
             zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
 
-            ! masks
+            ! Mask for ice presence (1) or absence (0)
             zmsk00x(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
             zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
 
-            ! switches 
+            ! Mask for lots of ice (1) or little ice (0)
             IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
             ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
             IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
-            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
+            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
 
          END DO
@@ -473,10 +459,8 @@
                
                ! delta* at T points (following Lemieux and Dupont, GMD 2020)
-               zdeltat_star = MAX( zdeltat, rn_delta_min )
-               
-               IF ( ln_smooth_vp ) zdelta_star = zdeltat + rn_delta_min
-               
+               zdeltastar_t(ji,jj) = zdeltat + rn_creepl
+              
                ! P/delta at T-points
-               zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltat_star
+               zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj)
                
                ! Temporary zzt and zet factors at T-points
@@ -486,5 +470,5 @@
             END DO
          END DO
-
+         
          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. )
 
@@ -502,5 +486,5 @@
          
          CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. )
-         
+
          !---------------------------------------------------
          ! -- Ocean-ice drag and Coriolis RHS contributions
@@ -799,5 +783,13 @@
                         !---------------
                         DO ji = 2, jpi - 1    
+                        
                            u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) )
+                           
+                           ! velocity masking for little-ice and no-ice cases
+                           ! put 0.01 of ocean current, appropriate choice to avoid too much spreading of the ice edge
+                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        & 
+                              &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             &
+                                          + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    )
+                           
                         END DO
 
@@ -871,4 +863,10 @@
                         DO jj = 2, jpj - 1
                            v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) )
+
+                           ! mask velocity for no-ice and little-ice cases                           
+                           v_ice(ji,jj) =   zmsk00y(ji,jj)                                        & 
+                              &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)             &
+                                          + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    )
+
                         END DO
 
@@ -897,5 +895,5 @@
                   DO jj = 2, jpj - 1
                      DO ji = 2, jpi - 1
-                        zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) ! * zmask15 ---> MV TEST mask at 15% concentration
+                        zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1)
                      END DO
                   END DO
@@ -975,5 +973,5 @@
       !------------------------------------------------------------------------------!
       !
-      ! OPT: subroutinize ?
+      ! MV OPT: subroutinize ?
             
       DO jj = 1, jpj - 1
@@ -997,4 +995,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  )  &
@@ -1013,5 +1013,7 @@
             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
+            
+            !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
+            pdelta_i(ji,jj) = zdelta + rn_creepl
 
          END DO
@@ -1019,15 +1021,5 @@
       
       CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
-      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
-      
-      ! --- Store the stress tensor for the next time step --- !
-      ! MV OPT: Are they computed at the end of the tucking fime step ???
-      ! MV Is stress at previous time step needed for VP, normally no, because they equation is not tucking fime dependent!!!
-      !
-      pstress1_i (:,:) = zs1 (:,:)
-      pstress2_i (:,:) = zs2 (:,:)
-      pstress12_i(:,:) = zs12(:,:)
-      !
-
+      
       !------------------------------------------------------------------------------!
       !
@@ -1035,13 +1027,57 @@
       !
       !------------------------------------------------------------------------------!
-      ! MV OPT: subroutinize
-      !
-      ! --- Ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
+      !
+      ! MV OPT: subroutinize ?
+      !
+      !----------------------------------
+      ! --- Recompute stresses if needed 
+      !----------------------------------
+      !
+      ! ---- Sea ice stresses at T-points
+      IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
+      
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
+               zp_deltastar_t(ji,jj)   =   strength(ji,jj) / pdelta_i(ji,jj) 
+               zfac                    =   zp_deltastar_t(ji,jj) 
+               zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
+               zs2(ji,jj)              =   zfac * z1_ecc2 * zten_i(ji,jj)
+               zs12(ji,jj)             =   zfac * z1_ecc2 * pshear_i(ji,jj)
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
+      
+      ENDIF
+      
+      ! ---- s12 at F-points      
+      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
+
+         DO jj = 1, jpj - 1
+            DO ji = 1, jpi - 1
+            
+               ! P/delta* at F points
+               zp_deltastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) )
+               
+               ! s12 at F-points 
+               zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj)
+               
+            END DO
+         END DO
+
+         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
+         
+      ENDIF
+      !
+      !-----------------------
+      ! --- Store diagnostics
+      !-----------------------
+      !
+      ! --- Ice-ocean, ice-atm. & ice-ocean bottom (landfast) stresses --- !
       IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
          & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
-
          !
          CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
-            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
          !
          CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
@@ -1059,38 +1095,70 @@
       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
+         !
+         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags.
+         ! Definitions following Coon (1974) and Feltham (2008)
+         !
+         ! sigma1, sigma2, sigma12 are useful (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
+         ! however these are NOT stress tensor components, neither stress invariants, nor stress principal components
+         !
+         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !         
          DO jj = 2, jpj - 1
             DO ji = 2, jpi - 1
-            
-               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) )          ! 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 )
-               
+               ! Stress invariants
+               zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp                                        ! 1st invariant, aka average normal stress aka negative pressure
+               zsig_II(ji,jj)   =   SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) )  ! 2nd invariant, aka maximum shear stress               
             END DO
          END DO
 
-         CALL lbc_lnk_multi( 'icedyn_rhg_vp', 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 )
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
+         
+         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 --- !
+      ! These are used to plot the normalized yield curve (Lemieux & Dupont, GMD 2020)
+      ! To plot the yield curve and evaluate physical convergence, they have two recommendations
+      ! Recommendation 1 : Use ice strength, not replacement pressure
+      ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765)
+      ! R2 means we need to recompute stresses
+      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
+         !
+         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )         
+         !         
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
+               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
+               !                        and **deformations** at current iterates
+               !                        following Lemieux & Dupont (2020)
+               zfac             =   zp_deltastar_t(ji,jj)
+               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(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), T-point
+               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                              ! 1st invariant
+               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd invariant
+
+               ! Normalized  principal stresses (used to display the ellipse)
+               z1_strength      =   1._wp / strength(ji,jj)
+               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  + zsig_II(ji,jj) ) * z1_strength
+               zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) - zsig_II(ji,jj) ) * z1_strength
+            END DO
+         END DO
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
+         !      
+         !
+         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
+         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
+
+         DEALLOCATE( zsig1_p , zsig2_p )
          
       ENDIF
@@ -1098,9 +1166,8 @@
       ! --- SIMIP, terms of tendency for momentum equation  --- !
       IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
-         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
-         !
-!!!!!!!!!         ATTENTION LA FORCE INTERNE DOIT ETTRE RECALCIULEEE ICCI !!!!!!!!!!!!!!!!
+         & iom_use('corstrx') .OR. iom_use('corstry') ) THEN
+         !
          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
-            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
+            &                                 zCorU, 'U', -1., zCorV, 'V', -1. )
 
          CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
@@ -1108,8 +1175,35 @@
          CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
          CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+      ENDIF
+      
+      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
+
+
+         ! Recalculate internal forces (divergence of stress tensor)
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
+               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e1e2u(ji,jj)
+               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e1e2v(ji,jj)
+            END DO
+         END DO
+            
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
+         
          CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
          CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
+         
       ENDIF
 
+      ! --- Ice & snow mass and ice area transports
       IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
          & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
@@ -1137,6 +1231,6 @@
 
          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
-            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
-            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+            &                                 zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
+            &                                 zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
 
          CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
@@ -1172,6 +1266,4 @@
                   &       prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV )
    
-!                  CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, &
-!                          &        zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE rhg_cvg_vp  ***
@@ -1302,60 +1394,4 @@
    
 
-
-   SUBROUTINE rhg_vp_rst( cdrw, kt )
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE rhg_vp_rst  ***
-      !!                     
-      !! ** Purpose :   Read or write RHG file in restart file
-      !!
-      !! ** Method  :   use of IOM library
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
-      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
-      !
-      INTEGER  ::   iter            ! local integer
-      INTEGER  ::   id1, id2, id3   ! local integers
-      !!----------------------------------------------------------------------
-      !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
-         !                                   ! ---------------
-         IF( ln_rstart ) THEN                   !* Read the restart file
-            !
-            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
-            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
-            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
-            !
-            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
-               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
-               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
-               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
-            ELSE                                     ! start rheology from rest
-               IF(lwp) WRITE(numout,*)
-               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
-               stress1_i (:,:) = 0._wp
-               stress2_i (:,:) = 0._wp
-               stress12_i(:,:) = 0._wp
-            ENDIF
-         ELSE                                   !* Start from rest
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
-            stress1_i (:,:) = 0._wp
-            stress2_i (:,:) = 0._wp
-            stress12_i(:,:) = 0._wp
-         ENDIF
-         !
-      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-         !                                   ! -------------------
-         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
-         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
-         !
-         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
-         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
-         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
-         !
-      ENDIF
-      !
-   END SUBROUTINE rhg_vp_rst
-
    
 #else
