Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90	(revision 13493)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90	(revision 13522)
@@ -150,10 +150,21 @@
    !
    !                                     !!** ice-rheology namelist (namdyn_rhg) **
+   ! -- evp
    LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
-   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
-   REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve
+   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit !!! common with vp
+   REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve !!! common with vp
    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 
+   ! -- vp
+   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)
+   REAL(wp), PUBLIC ::   rn_relaxv_vp     !: V-relaxation factor 
+   REAL(wp), PUBLIC ::   rn_uerr_max_vp   !: maximum velocity error, above which a forcing error is considered and solver is stopped
+   REAL(wp), PUBLIC ::   rn_uerr_min_vp   !: minimum velocity error, beyond which convergence is assumed
+   INTEGER , PUBLIC ::   nn_cvgchk_vp     !: Number of iterations every each convergence is checked
    !
    !                                     !!** ice-advection namelist (namdyn_adv) **
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 13493)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90	(revision 13522)
@@ -37,5 +37,5 @@
 
    ! ** namelist (namrhg) **
-   LOGICAL ::   ln_rhg_EVP       ! EVP rheology
+   LOGICAL ::   ln_rhg_EVP       ! EVP rheology ! MV: why is it declared here while all other parameters are declared in ice.F90
    LOGICAL ::   ln_rhg_VP        ! EVP rheology
    !
@@ -120,7 +120,8 @@
       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, ln_rhg_chkcvg, &
-    &                       ln_rhg_VP 
+      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
       !!-------------------------------------------------------------------
+      
       !
       REWIND( numnam_ice_ref )         ! Namelist namdyn_rhg in reference namelist : Ice dynamics
@@ -139,11 +140,21 @@
          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,*) '         creep limit                                       rn_creepl     = ', rn_creepl     ! common with VP
+         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc        ! common with VP
          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
+         WRITE(numout,*) '      check convergence of rheology                        ln_rhg_chkcvg = ', ln_rhg_chkcvg ! common with VP
          WRITE(numout,*) ''    
-         WRITE(numout,*) '      rheology VP     (icedyn_rhg_lsr)                     ln_rhg_VP     = ', ln_rhg_VP
+         WRITE(numout,*) '      rheology VP     (icedyn_rhg_VP)                      ln_rhg_VP     = ', ln_rhg_VP
+         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
+         WRITE(numout,*) '         relaxation factor for v                           rn_relaxv_vp  = ', rn_relaxv_vp
+         WRITE(numout,*) '         maximum error on velocity                         rn_uerr_max_vp = ', rn_uerr_max_vp
+         WRITE(numout,*) '         velocity to decide convergence                    rn_uerr_min_vp = ', rn_uerr_min_vp
+         WRITE(numout,*) '         iteration step for convergence check              nn_cvgchk_vp   = ', nn_cvgchk_vp
+                           
       ENDIF
       !
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 13493)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13522)
@@ -1,33 +1,39 @@
 ! TODOLIST
-
-!
-! Code residual norm of the linear system
-! Verify convergence check
 !
 ! Define all symbols
-! Check uses
+! - 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
+! 
+!
+! Commit
 ! Compile
 !
-! Check
-! --- lbclnks --> what is the policy ?
-! --- boundary conditions --> how to enforce them - is the fmask strategy sufficient ?
-! --- change of independent term in the U-V solvers ?
-! --- loop ji then jj in the V-solver, is it well-written
+! 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
 !
-! Subroutinize tridiagonal systems ?
+! 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
+! - Tensile strength
+! - agrif
+! - bdy
+!
 ! Write documentation
 !
-
-! Convergence quality criteria
-! - average KE (for outer loop iterations) (see Lemieux et al, 2009)
-! - residual of the linear system (for inner loop iterations)
-! - velocity difference ?
-
-
-! add bottom drag
-! add tensile strength
-
-! quality control - compare to mitGCM
-! quality control - comparing EVP and VP
+! Optimize
+! - subroutinize common parts (diagnostics)
+! - namelist: rationalize common parameters EVP/VP
+
 
 MODULE icedyn_rhg_vp
@@ -75,10 +81,13 @@
    PUBLIC   rhg_vp_rst       ! called by icedyn_rhg.F90
 
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-
    !! for convergence tests
    INTEGER ::   ncvgid   ! netcdf file id
-   INTEGER ::   nvarid   ! netcdf variable id
+   INTEGER ::   nvarid_ucvg   ! netcdf variable id
+   INTEGER ::   nvarid_ures
+   INTEGER ::   nvarid_vres
+   INTEGER ::   nvarid_velres
+   INTEGER ::   nvarid_udif
+   INTEGER ::   nvarid_vdif
+   INTEGER ::   nvarid_mke
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
    !!----------------------------------------------------------------------
@@ -87,4 +96,5 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
+   
 CONTAINS
 
@@ -95,24 +105,29 @@
       !!                             VP-LSR-C-grid
       !!
-      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
+      !! ** Purpose : determines sea ice drift from wind stress, ice-ocean
       !!  stress and sea-surface slope. Internal forces assume viscous-plastic rheology (Hibler, 1979)
+      !! 
+      !! ** Method
       !!  
       !!  The resolution algorithm  follows from Zhang and Hibler (1997) and Losch (2010)
+      !!  with elements from Lemieux and Tremblay (2008) and Lemieux and Tremblay (2009)
       !!  
-      !!  It is based on a linear slope relaxation algorithm (LSR)
-      !!
-      !!  * Time stepping
-      !!
-      !!  The idea is to arrange the momentum equations as follows:
+      !!  The components of the momentum equations are arranged following the ideas of Zhang and Hibler (1997)
+      !!
       !!  f1(u) = g1(v)
       !!  f2(v) = g2(v)
       !!
-      !!  and use two sub-time-step levels (outer and inner) to solve them.
-      !!
-      !!  The right-hand side (RHS, g1 & g2) are expressed explicitly, and the left-hand side (LHS, f1, f2) implicitly
-      !!
-      !!  The inner loop makes ~1500 steps max to solve a linear system problem, targetting to solve the implicit terms
-      !!  The outer loop makes ~10 steps and updates the explicit terms
-      !!  Outer loop uses a modified euler time step, with uC = 0.5 * ( uC + u )
+      !!  The right-hand side (RHS) is explicit
+      !!  The left-hand side (LHS) is implicit
+      !!  Coriolis is part of explicit terms, whereas ice-ocean drag is implicit
+      !!
+      !!  Two iteration levels (outer and inner loops) are used to solve the equations
+      !!
+      !!  The outer loop (OL, typically 10 iterations) is there to deal with the (strong) non-linearities in the equation
+      !!
+      !!  The inner loop (IL, typically 1500 iterations) is there to solve the linear problem with a line-successive-relaxation algorithm
+      !!
+      !!  The velocity used in the non-linear terms uses a "modified euler time step" (not sure its the correct term), 
+      !!! with uk = ( uk-1 + uk-2 ) / 2.
       !!  
       !!  * Spatial discretization 
@@ -143,86 +158,80 @@
       !! ** Notes   : 
       !!             
-      !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010. 
+      !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010., Lemieux et al., 2008, 2009, ...  
       !!             
       !!             
       !!-------------------------------------------------------------------
       !!
-!!!      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      !
-!!!      !!
-!!!      INTEGER ::   ji, jj       ! dummy loop indices
-!!!      INTEGER ::   jter         ! local integers
-!!!      !
-!!!      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
-!!!      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
-!!!      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) ::   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) ::   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) ::   zbeta                           ! beta coef from Kimmritz 2017
-!!!      !
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   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) ::   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
-!!!      !                                                                 !    ice bottom surface if ice is embedded   
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
-!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
+      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      !
+      !!
+      LOGICAL ::   ll_u_iterate, ll_viterate   ! continue iteration or not
+      !
+      INTEGER ::   ji, jj, jn          ! dummy loop indices
+      INTEGER ::   jter, i_out, i_inn  ! 
+      INTEGER ::   ji_min, jj_min      !
+      INTEGER ::   nn_zebra_vp         ! number of zebra steps
+      INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)      
+      !
+      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
+      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
+      REAL(wp) ::   zglob_area                                          ! global ice area for diagnostics
+      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) ::   zp_deltastar_f
+      REAL(wp) ::   zfac, zfac1, zfac2, zfac3
+      REAL(wp) ::   zt12U, zt11U, zt22U, zt21U, zt122U, zt121U
+      REAL(wp) ::   zt12V, zt11V, zt22V, zt21V, zt122V, zt121V
+      REAL(wp) ::   zAA3, zw, zuerr_max, zverr_max
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
+      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! Acceleration term contribution to RHS
+      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) ::   zzt, zet                        ! Viscosity pre-factors at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zef                             ! Viscosity pre-factor at F point
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmt                             ! Mass per unit area at t-point
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! Coriolis factor (m*f) at t-point 
+      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) ::   zu_c, zv_c                      ! "current" ice velocity (m/s), average of previous two OL iterates
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_cV, zv_cU                    ! "current" u (v) ice velocity at V (U) point (m/s)
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_b, zv_b                      ! velocity at previous sub-iterate
+      REAL(wp), DIMENSION(jpi,jpj) ::   zuerr, zverr                    ! absolute U/Vvelocity difference between current and previous sub-iterates
+      !
+      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
+      !                                                                 !    ice bottom surface if ice is embedded   
+      REAL(wp), DIMENSION(jpi,jpj) ::   zCwU, zCwV                      ! ice-ocean drag pre-factor (rho*c*module(u))
+      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi_rhsu, ztauy_oi_rhsv    ! ice-ocean stress RHS contribution at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsu, zs2_rhsu, zs12_rhsu   ! internal stress contributions to RHSU
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsv, zs2_rhsv, zs12_rhsv   ! internal stress contributions to RHSV
+      REAL(wp), DIMENSION(jpi,jpj) ::   zf_rhsu, zf_rhsv                ! U- and V- components of internal force RHS contributions
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrhsu, zrhsv                    ! U and V RHS 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zAU, zBU, zCU, zDU, zEU         ! Linear system coefficients, U equation
+      REAL(wp), DIMENSION(jpi,jpj) ::   zAV, zBV, zCV, zDV, zEV, zFV    ! Linear system coefficients, V equation
+      REAL(wp), DIMENSION(jpi,jpj) ::   zFU, zFU_prime, zBU_prime       ! Rearranged linear system coefficients, U equation
+      REAL(wp), DIMENSION(jpi,jpj) ::   zFV, zFV_prime, zBV_prime       ! Rearranged linear system coefficients, V equation
 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
 !!!      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), 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
-!!!
-!!!      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
-
-         !! --- check convergence
-         
-         
-         !! --- LSR arrays
-         REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                           ! mask at F points for the ice
-
-         REAL(wp), DIMENSION(jpi,jpj) ::   zaU      , zaV                   ! ice fraction on U/V points
-         REAL(wp), DIMENSION(jpi,jpj) ::   zmassU   , zmassV                ! ice mass per unit area at U/V points (kg/m2)
-         
-         REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice_C , zv_ice_C              ! ice velocities averaged between before and mid outer sub-iteration
-         REAL(wp), DIMENSION(jpi,jpj) ::   uTmp    , vTmp                   ! uTmp  = velocity at previous *inner* interation 
-
-         REAL(wp), DIMENSION(jpi,jpj) ::   z_etat  , z_zetat                ! VP viscosities at T points
-         REAL(wp), DIMENSION(jpi,jpj) ::   z_etaf                           ! VP eta viscosity at F point
-         
-         REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsu, zs2_rhsu, zs12_rhsu    ! Internal stress contribution to RHS of U-equation
-         REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsv, zs2_rhsv, zs12_rhsv    ! Internal stress contribution to RHS of V-equation
-         REAL(wp), DIMENSION(jpi,jpj) ::   zf_rhsu, zf_rhsv                 ! U- and V- components of internal forces
-         
-!!!      !! --- SIMIP diags
+     !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
+      !
+      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
+      !! --- 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)
@@ -231,34 +240,28 @@
       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)      
-                  
-      !! --- Convergence diags
-      REAL(wp), ALLOCATABLE, DIMENSION(:) ::  zdiag_meanke(:)
-      REAL(wp), ALLOCATABLE, DIMENSION(:) ::  zdiag_resnor(:)
-      REAL(wp), ALLOCATABLE, DIMENSION(:) ::  zdiag_veldif(:)
-      
+                        
       !!----------------------------------------------------------------------------------------------------------------------
 
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)'
-      
-      ! Namelist parameters to outsource (once we all are super happy)
-      nn_nouter_vp      = 3            ! Number of outer iterations (mitGCM 10) 
-      nn_ninner_vp      = 1500         ! Number of inner-loop iterations (mitGCM 1500)
-      rn_deltamin       = 2.0e-9       ! minimum delta value following Lemieux and Dupont, GMD 2020                          
-      ln_visc_VP_smooth = .FALSE.      ! zeta viscosity smoothing (if yes, tanh function of Lemieux and Tremblay JGR 2009)
-
-      ! Linear system parameters
-      ln_zebra          = .FALSE.      ! do two passages for the linear system problem, one every odd j-band, one every even one (recommended by Martin Losch)
-      IF ( ln_zebra ) THEN; nn_zebra_step = 2; ELSE; nn_zebra_step = 1; ENDIF 
-
-      rn_relaxfacU = 0.95 ! U-relaxation factor for linear problem (Losch 0.95, Zhang 1.8)
-      rn_relaxfacV = 0.95 ! V-relaxation factor   
-      
-      rn_uerr_max = 0.8_wp ! maximum error on velocity, above which a forcing error is considered
-      rn_uerr_min = 1.0d-4 ! velocity difference beyond which the linear system procedure is stopped
-      nn_lsr_cvgcheck = 5  ! iterations every each convergence is checked
-
+            
       !------------------------------------------------------------------------------!
-      ! 1) Initialization
+      !
+      ! --- 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
+      
+      IF ( ln_zebra_vp ) THEN; nn_nzebra_vp = 2; ELSE; nn_nzebra_vp = 1; ENDIF 
+      
+      nn_nvp = nn_out_vp * nn_inn_vp ! maximum number of iterations
+      
       zrhoco = rau0 * rn_cio 
 
@@ -272,24 +275,14 @@
       zs12(:,:) = pstress12_i(:,:)
       
-      ! Convergence checks
-      IF( ln_rhg_chkcvg ) THEN 
-      
+      ! Initialise convergence checks
+      IF( ln_rhg_chkcvg ) THEN
+      
+         ! ice area for global mean kinetic energy
          zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) ) ! global ice area (km2)
-         z1_32      = 1._wp / 32._wp
-
-         ! --- mean kinetic energy      
-         ALLOCATE( zdiag_meanke(nn_outer_vp*nn_inner_vp) )  ! mean kinetic energy          
-         ALLOCATE( zdiag_resnor(nn_outer_vp*nn_inner_vp) )  ! residual norm of the linear system  
-         ALLOCATE( zdiag_u_idif(nn_outer_vp*nn_inner_vp) )  ! velocity difference between two iterations
-         ALLOCATE( zdiag_v_idif(nn_outer_vp*nn_inner_vp) )  ! velocity difference between two iterations
-         
-         zdiag_resnor(:) = 0._wp
-         zdiag_meanke(:) = 0._wp
-         zdiag_u_idif(:) = 0._wp
-         zdiag_v_idif(:) = 0._wp
          
       ENDIF
 
-      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
+      ! Landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
+      ! MV: Not working yet...
       IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
       ELSE                         ;   zkt = 0._wp
@@ -297,19 +290,21 @@
 
       !------------------------------------------------------------------------------!
-      ! 2) Time-independent quantities
+      !
+      ! --- Time-independent quantities
+      !
       !------------------------------------------------------------------------------!
       
       CALL ice_strength ! strength at T points
       
-      !--------
-      ! F-mask      
-      !--------
+      !------------------------------
+      ! -- F-mask       (code from EVP)
+      !------------------------------
       ! MartinV: 
       ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip)
-      ! I don't know exactly where it should be applied in here
-      
-      ! ocean/land mask
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1      ! NO vector opt.
+      ! I am not sure the same recipe applies here
+      
+      ! - ocean/land mask
+      DO jj = 1, jpj - 1
+         DO ji = 1, jpi - 1
             zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
          END DO
@@ -318,6 +313,7 @@
 
       ! Lateral boundary conditions on velocity (modify zfmask)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      ! Can be computed once for all, at first time step, for all rheologies
+      DO jj = 2, jpj - 1
+         DO ji = 2, jpi - 1   ! vector opt.
             IF( zfmask(ji,jj) == 0._wp ) THEN
                zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
@@ -326,18 +322,18 @@
          END DO
       END DO
-      DO jj = 2, jpjm1
+      DO jj = 2, jpj - 1
          IF( zfmask(1,jj) == 0._wp ) THEN
             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( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) )
+            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpi - 1,jj,1), umask(jpi,jj-1,1) ) )
         ENDIF
       END DO
-      DO ji = 2, jpim1
+      DO ji = 2, jpi - 1
          IF( zfmask(ji,1) == 0._wp ) THEN
             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( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) )
+            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpj - 1,1) ) )
          ENDIF
       END DO
@@ -345,12 +341,16 @@
       CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp )
       
-      !---------------------------------------------------------------------------------------------------------
-      ! 3) Time-independent contributions of acceleration, ocean drag, coriolis, atmospheric drag, surface tilt
-      !---------------------------------------------------------------------------------------------------------
-      ! 
-      ! --- Compute all terms & factors independent of velocities, or only depending on velocities at previous time step
-      
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
+      !----------------------------------------------------------------------------------------------------------
+      ! -- Time-independent pre-factors for acceleration, ocean drag, coriolis, atmospheric drag, surface tilt
+      !----------------------------------------------------------------------------------------------------------
+      ! Compute all terms & factors independent of velocities, or only depending on velocities at previous time step      
+      
+      ! sea surface height
+      !    embedded sea ice: compute representative ice top surface
+      !    non-embedded sea ice: use ocean surface for slope calculation
+      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
+      
+      DO jj = 2, jpj - 1
+         DO ji = 2, jpi - 1
 
             ! Ice fraction at U-V points
@@ -359,12 +359,12 @@
 
             ! Snow and ice mass at U-V points
-            zm1             = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
+            zmt(ji,jj)      = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
             zm2             = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
             zm3             = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
             
-            zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
-            zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+            zmassU          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zmassV          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
                           
-            ! Mass per unit time
+            ! Mass per unit area divided by time step
             zmassU_t(ji,jj) = zmassU * r1_rdtice
             zmassV_t(ji,jj) = zmassV * r1_rdtice
@@ -379,5 +379,5 @@
 
             ! Coriolis factor at T points (m*f)
-            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
+            zmf(ji,jj)      = zmt(ji,jj) * ff_t(ji,jj)
             
             ! Wind stress
@@ -393,5 +393,5 @@
             zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
 
-            ! switches
+            ! switches 
             IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
             ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
@@ -403,54 +403,43 @@
             
       !------------------------------------------------------------------------------!
-      ! 4) Start outer loop, update velocities
+      !
+      ! --- Start outer loop
+      !
       !------------------------------------------------------------------------------!
       
-      zu_ice_C(:,:) = u_ice(:,:)
-      zv_ice_C(:,:) = v_ice(:,:)
-      
-      i_iter = 0
-
-      DO i_out = 1, nn_nouter_vp
-      
-         i_iter = i_iter + 1
-      
-         ! EVP CODE
-         ! l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
-         !
-         ! 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
-         
-         ! Zhang and Hibler 97 time stepping
-         ! identified as the most rapid test with mitGCM and recommendation by Martin Losch
-         zu_ice_C(:,:) = 0.5_wp * ( u_ice(:,:) + zu_ice_C(:,:) )
-         zv_ice_C(:,:) = 0.5_wp * ( v_ice(:,:) + zv_ice_C(:,:) )
-         
-      !------------------------------------------------------------------------------!
-      ! 5) Right-hand side of linear problem
-      !------------------------------------------------------------------------------!
-
-         !      
+      zu_c(:,:) = u_ice(:,:)
+      zv_c(:,:) = v_ice(:,:)
+      
+      jter = 0
+
+      DO i_out = 1, nn_nout_vp
+               
+         ! Velocities used in the non linear terms are the average of the past two iterates
+         ! u_it = 0.5 * ( u_{it-1} + u_{it-2})
+         ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009)
+         zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) )
+         zv_c(:,:) = 0.5_wp * ( v_ice(:,:) + zv_c(:,:) )
+         
+         !------------------------------------------------------------------------------!
+         !
+         ! --- Right-hand side (RHS) of the linear problem
+         !
+         !------------------------------------------------------------------------------!
          ! In the outer loop, one needs to update all RHS terms
          ! with explicit velocity dependencies (viscosities, coriolis, ocean stress)
-         !
-      
-         !--------------------------------------------------------
-         ! 5.1 Strain rates, viscosities and replacement pressure
-         !--------------------------------------------------------
+         ! as a function of uc
+         !
+      
+         !------------------------------------------
+         ! -- Strain rates, viscosities and P/Delta
+         !------------------------------------------
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
-            DO ji = 1, jpim1
+         DO jj = 1, jpj - 1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+            DO ji = 1, jpi - 1
 
                ! shear at F points
-               zds(ji,jj) = ( ( zu_ice_C(ji,jj+1) * r1_e1u(ji,jj+1) - zu_ice_C(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
-                  &         + ( zv_ice_C(ji+1,jj) * r1_e2v(ji+1,jj) - zv_ice_C(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               zds(ji,jj) = ( ( zu_c(ji,jj+1) * r1_e1u(ji,jj+1) - zu_c(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+                  &         + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
                   &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
 
@@ -458,8 +447,8 @@
          END DO
          
-         CALL lbc_lnk( 'icedyn_rhg_vp', 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
+         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan
+
+         DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
+            DO ji = 2, jpi - 1 ! 
 
                ! shear**2 at T points (doc eq. A16)
@@ -469,12 +458,12 @@
               
                ! divergence at T points
-               zdiv  = ( e2u(ji,jj) * zu_ice_C(ji,jj) - e2u(ji-1,jj) * zu_ice_C(ji-1,jj)   &
-                  &    + e1v(ji,jj) * zv_ice_C(ji,jj) - e1v(ji,jj-1) * zv_ice_C(ji,jj-1)   &
+               zdiv  = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj)   &
+                  &    + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1)   &
                   &    ) * r1_e1e2t(ji,jj)
                zdiv2 = zdiv * zdiv
                
                ! tension at T points
-               zdt  = ( ( zu_ice_C(ji,jj) * r1_e2u(ji,jj) - zu_ice_C(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-                  &   - ( zv_ice_C(ji,jj) * r1_e1v(ji,jj) - zv_ice_C(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               zdt  = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &   - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
                   &   ) * r1_e1e2t(ji,jj)
                zdt2 = zdt * zdt
@@ -486,7 +475,5 @@
                zdeltat_star = MAX( zdeltat, rn_delta_min )
                
-!               IF ( ln_visc_VP_smooth ) THEN ! to code, other regularization from Lemieux and Tremblay (2009)
-!                   zdeltat_star = ...
-!               ENDIF
+               IF ( ln_smooth_vp ) zdelta_star = zdeltat + rn_delta_min
                
                ! P/delta at T-points
@@ -494,5 +481,5 @@
                
                ! Temporary zzt and zet factors at T-points
-               zzt(ji,jj)     = zp_deltastar_t * r1_e1e2t(ji,jj)
+               zzt(ji,jj)     = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)
                zet(ji,jj)     = zzt(ji,jj)     * z1_ecc2 
                           
@@ -500,10 +487,8 @@
          END DO
 
-         CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. )
-         CALL lbc_lnk( 'icedyn_rhg_vp', zzt , 'T', 1. )
-         CALL lbc_lnk( 'icedyn_rhg_vp', zet, 'T', 1. )
-
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. )
+
+         DO jj = 1, jpj - 1
+            DO ji = 1, jpi - 1
                
                ! P/delta* at F points
@@ -519,19 +504,19 @@
          
          !---------------------------------------------------
-         ! 5.2 Ocean-ice drag and Coriolis RHS contributions
+         ! -- Ocean-ice drag and Coriolis RHS contributions
          !---------------------------------------------------
          
-         DO jj = 2, jpjm1
-             DO ji = fs_2, fs_jpim1
+         DO jj = 2, jpj - 1
+             DO ji = 2, jpi - 1
                 
                 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
-                zu_ice_CV(ji,jj)     = 0.25_wp * ( zu_ice_C(ji,jj) + zu_ice_C(ji-1,jj) + zu_ice_C(ji,jj+1) + zu_ice_C(ji-1,jj+1) ) * vmask(ji,jj,1)
-                zv_ice_CU(ji,jj)     = 0.25_wp * ( zv_ice_C(ji,jj) + zv_ice_C(ji,jj-1) + zv_ice_C(ji+1,jj) + zv_ice_C(ji+1,jj-1) ) * umask(ji,jj,1)
-
+                zu_cV            = 0.25_wp * ( zu_c(ji,jj) + zu_c(ji-1,jj) + zu_c(ji,jj+1) + zu_c(ji-1,jj+1) ) * vmask(ji,jj,1)
+                zv_cU            = 0.25_wp * ( zv_c(ji,jj) + zv_c(ji,jj-1) + zv_c(ji+1,jj) + zv_c(ji+1,jj-1) ) * umask(ji,jj,1)
+                
                 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
-                zCwU(ji,jj)          = zaU(ji,jj) * zrhoco * SQRT( ( zu_ice_C (ji,jj) - u_oce (ji,jj) ) * ( zu_ice_C (ji,jj) - u_oce (ji,jj) )  &
-                  &                                              + ( zv_ice_CU(ji,jj) - v_oceU(ji,jj) ) * ( zv_ice_CU(ji,jj) - v_oceU(ji,jj) ) )
-                zCwV(ji,jj)          = zaV(ji,jj) * zrhoco * SQRT( ( zv_ice_C (ji,jj) - v_oce (ji,jj) ) * ( zv_ice_C (ji,jj) - v_oce (ji,jj) )  &
-                  &                                              + ( zu_ice_CV(ji,jj) - u_oceV(ji,jj) ) * ( zu_ice_CV(ji,jj) - u_oceV(ji,jj) ) )
+                zCwU(ji,jj)          = zaU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  &
+                  &                                              + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
+                zCwV(ji,jj)          = zaV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  &
+                  &                                              + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
                  
                 !--- Ocean-ice drag contributions to RHS 
@@ -542,63 +527,66 @@
                 ! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done
                 zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  &
-                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_ice_C(ji  ,jj) + e1v(ji  ,jj-1) * zv_ice_C(ji  ,jj-1) )  &
-                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_ice_C(ji+1,jj) + e1v(ji+1,jj-1) * zv_ice_C(ji+1,jj-1) ) )
+                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_c(ji  ,jj) + e1v(ji  ,jj-1) * zv_c(ji  ,jj-1) )  &
+                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) )
                            
                 ! --- V-component of Coriolis Force (energy conserving formulation)
                 zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  &
-                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_ice_C(ji,jj  ) + e2u(ji-1,jj  ) * zu_ice_C(ji-1,jj  ) )  &
-                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_ice_C(ji,jj+1) + e2u(ji-1,jj+1) * zu_ice_C(ji-1,jj+1) ) )
+                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_c(ji,jj  ) + e2u(ji-1,jj  ) * zu_c(ji-1,jj  ) )  &
+                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) )
          
              END DO
          END DO
          
-         ! ??? lbclnk ???
-
-         !--------------------------------------
-         ! 5.3 Internal stress RHS contribution
-         !--------------------------------------
-
-         ! --- Stress contributions at t-points         
+         ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system, 
+         ! so there is no need for lbclnk on drag and coriolis
+
+         !-------------------------------------
+         ! -- Internal stress RHS contribution
+         !-------------------------------------
+
+         ! --- Stress contributions at T-points         
          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 ji = 2, jpi ! 
             
                ! sig1 contribution to RHS of U-equation at T-points
-               zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_ice_C(ji,jj) - e1v(ji,jj-1)    * zv_ice_C(ji,jj-1) - 1.0_wp )
+               zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_c(ji,jj) - e1v(ji,jj-1)    * zv_c(ji,jj-1) - 1.0_wp )
                                             
                ! sig2 contribution to RHS of U-equation at T-points            
-               zs2_rhsu(ji,jj) = - zet(ji,jj) * ( r1_e1v(ji,jj) * zv_ice_C(ji,jj) - r1_e1v(ji,jj-1) * zv_ice_C(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) 
+               zs2_rhsu(ji,jj) = - zet(ji,jj) * ( r1_e1v(ji,jj) * zv_c(ji,jj) - r1_e1v(ji,jj-1) * zv_c(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) 
 
                ! sig1 contribution to RHS of V-equation at T-points
-               zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj) * zu_ice_C(ji,jj)    - e2u(ji-1,jj)    * zu_ice_C(ji-1,jj) - 1.0_wp )
+               zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj)    * zu_c(ji,jj) - e2u(ji-1,jj)    * zu_c(ji-1,jj) - 1.0_wp )
 
                ! sig2 contribution to RHS of V-equation  at T-points
-               zs2_rhsv(ji,jj) =   zet(ji,jj) * ( r1_e2u(ji,jj) * zu_ice_C(ji,jj) - r1_e2u(ji-1,jj) * zu_ice_C(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)
+               zs2_rhsv(ji,jj) =   zet(ji,jj) * ( r1_e2u(ji,jj) * zu_c(ji,jj) - r1_e2u(ji-1,jj) * zu_c(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)
 
             END DO
          END DO
          
-         ! ??? lbclnk ???
+         ! a priori, no lbclnk, because rhsu is only used in the inner domain
          
          ! --- Stress contributions at f-points         
-         ! Note from MartinV: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing
+         ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing
          ! My guess is that this is the way to enforce boundary conditions on strain rate tensor
 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
+         DO jj = 1, jpj - 1
+            DO ji = 1, jpi - 1
                
                ! sig12 contribution to RHS of U equation at F-points 
-               zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_ice_C(ji+1,jj) - r1_e2v(ji,jj) * zv_ice_C(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj)
+               zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) - r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj)
                
                ! sig12 contribution to RHS of V equation at F-points
-               zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_ice_C(ji,jj+1) - r1_e1u(ji,jj) * zu_ice_C(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj)
+               zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) - r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj)
 
             END DO
          END DO
          
-         ! ??? lbcblnk
-
-         ! --- Internal forces contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1               
+         ! a priori, no lbclnk, because rhsu are only used in the inner domain
+
+         ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002)
+         ! OPT: merge with next loop and use intermediate scalars for zf_rhsu
+         
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1               
                ! --- U component of internal force contribution to RHS at U points
                zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & 
@@ -606,4 +594,5 @@
                   &      +         r1_e2u(ji,jj) * ( e2t(ji+1,jj) * e2t(ji+1,jj) * zs2_rhsu(ji+1,jj) - e2t(ji,jj)   * e2t(ji,jj)   * zs2_rhsu(ji,jj) )     &
                   &      + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj)   * e1f(ji,jj)   * zs12_rhsu(ji,jj)  - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) ) 
+                  
                ! --- V component of internal force contribution to RHS at V points
                zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * &
@@ -611,14 +600,15 @@
                   &      +         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj)     * zs2_rhsv(ji,jj) )     &
                   &      + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj)   * e2f(ji,jj)   * zs12_rhsv(ji,jj)  - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) )
+                  
             END DO
          END DO
          
-         ! ??? lbclnk ???
-
          !---------------------------
-         ! 5.4 Sum RHS contributions
+         ! -- Sum RHS contributions
          !---------------------------
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
+         !
+         ! OPT: could use intermediate scalars to reduce memory access
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
             
                ! still miss ice ocean stress and acceleration contribution
@@ -629,9 +619,11 @@
          END DO
          
-      ! ??? lbclnk ???
+         ! inner domain calculations -> no lbclnk
      
-      !---------------------------------------------------------------------------------------!
-      ! 6) Linear system matrix
-      !---------------------------------------------------------------------------------------!
+         !---------------------------------------------------------------------------------------!
+         !
+         ! --- Linear system matrix
+         !
+         !---------------------------------------------------------------------------------------!
       
          ! Linear system matrix contains all implicit contributions 
@@ -641,135 +633,220 @@
          ! AU * u_{i-1,j} + BU * u_{i,j}   + CU * u_{i+1,j}
          !                = DU * u_{i,j-1} + EU * u_{i,j+1} + RHS 			(! my convention, not the same as ZH97 )         
-         !
-         !--------------------------------------
-         ! 6.1 Internal forces LHS contribution
-         !--------------------------------------
-         !
-         ! Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ?
-         ! Note 2: "T" factor calculations can be optimized by putting things out of the loop 
+         
+         ! MV Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ?
+         ! MV Note 2: "T" factor calculations can be optimized by putting things out of the loop 
          !         only zzt and zet are iteration-dependent, other only depend on scale factors
-         !
-         ! --- U-component
-         !
-         ! "T" factors (intermediate results)
-         !
-         zfac       = 0.5_wp * r1_e1e2u(ji,jj)
-         zfac1      =         zfac * e2u(ji,jj)
-         zfac2      =         zfac * r1_e2u(ji,jj)
-         zfac3      = 2._wp * zfac * r1_e1u(ji,jj)
-
-         zt12U      = - zfac1 * zzt(ji+1,jj)
-         zt11U      =   zfac1 * zzt(ji,jj)
-         
-         zt22U      = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj)
-         zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)
-         
-         zt122U     = - zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)
-         zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)
-
-         !
-         ! Linear system coefficients
-         !
-         zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj)
-         zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj)
-         zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj)
-         
-         zDU(ji,jj) =   zt121U * r1_e1u(ji,jj-1)
-         zEU(ji,jj) =   zt122U * r1_e1u(ji,jj+1)
-         
-         !
-         ! --- V-component
-         !
-         ! "T" factors (intermediate results)
-         !
-         zfac       = 0.5_wp * r1_e1e2v(ji,jj)
-         zfac1      =         zfac * e2v(ji,jj)
-         zfac2      =         zfac * r1_e1v(ji,jj)
-         zfac3      = 2._wp * zfac * r1_e2v(ji,jj)
-         
-         zt12V      = - zfac1 * zzt(ji,jj+1)
-         zt11V      =   zfac1 * zzt(ji,jj)
-         
-         zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1)
-         zt21V      = - zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)
-         
-         zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)
-         zt121V     = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)
-         
-         !
-         ! Linear system coefficients
-         !
-         zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1)
-         zBV(ji,jj) =  ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj)
-         zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1)
-         
-         zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here
-         zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj)
-                  
-         !---------------------------------
-         ! Ocean-ice drag LHS contribution
-         !---------------------------------
-         zBU(i,j) = zBU(i,j) + zCwU(ji,jj)
-
-         !-------------------------------
-         ! Acceleration LHS contribution
-         !-------------------------------
-         zBU(i,j) = zBU(i,j) + zmassU_t(ji,jj) 
+                  
+         DO ji = 2, jpi - 1 ! internal domain do loop
+            DO jj = 2, jpj - 1
+
+               !-------------------------------------
+               ! -- Internal forces LHS contribution
+               !-------------------------------------
+               !
+               ! --- U-component
+               !
+               ! "T" factors (intermediate results)
+               !
+               zfac       = 0.5_wp * r1_e1e2u(ji,jj)
+               zfac1      =         zfac * e2u(ji,jj)
+               zfac2      =         zfac * r1_e2u(ji,jj)
+               zfac3      = 2._wp * zfac * r1_e1u(ji,jj)
+
+               zt12U      = - zfac1 * zzt(ji+1,jj)
+               zt11U      =   zfac1 * zzt(ji,jj)
+         
+               zt22U      = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj)
+               zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)
+         
+               zt122U     = - zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)
+               zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)
+               
+               !
+               ! Linear system coefficients
+               !
+               zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj)
+               zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj)
+               zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj)
+         
+               zDU(ji,jj) =   zt121U * r1_e1u(ji,jj-1)
+               zEU(ji,jj) =   zt122U * r1_e1u(ji,jj+1)
+              
+               !
+               ! --- V-component
+               !
+               ! "T" factors (intermediate results)
+               !
+               zfac       = 0.5_wp * r1_e1e2v(ji,jj)
+               zfac1      =         zfac * e2v(ji,jj)
+               zfac2      =         zfac * r1_e1v(ji,jj)
+               zfac3      = 2._wp * zfac * r1_e2v(ji,jj)
+         
+               zt12V      = - zfac1 * zzt(ji,jj+1)
+               zt11V      =   zfac1 * zzt(ji,jj)
+         
+               zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1)
+               zt21V      = - zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)
+         
+               zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)
+               zt121V     = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)
+         
+               !
+               ! Linear system coefficients
+               !
+               zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1)
+               zBV(ji,jj) =  ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj)
+               zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1)
+         
+               zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here
+               zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj)
+                  
+               !-----------------------------------------------------
+               ! -- Ocean-ice drag and acceleration LHS contribution
+               !-----------------------------------------------------
+               zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj)
+               zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj)
+         
+            END DO
+         END DO
                
       !------------------------------------------------------------------------------!
-      ! 6) inner loop: solve linear system, check convergence
+      !
+      ! --- Inner loop: solve linear system, check convergence
+      !
       !------------------------------------------------------------------------------!
                
          ! Inner loop solves the linear problem .. requires 1500 iterations
-         doIterate4u = .TRUE.
-         doIterate4v = .TRUE.
-
-         DO i_inner = 1, nn_ninner_vp ! inner loop iterations
-         
-            ! initialize residual ...
-            
-            ...
-
-            IF ( doIterate4u .OR. doIterate4v ) THEN
-
-               zu_p(:,:) = u_ice(:,:) ! previous sub-iteration value
-               zv_p(:,:) = v_ice(:,:)
-
-                                          ! ---------------------------- !
-               IF ( doIterate4u ) THEN    ! --- Solve for u-velocity --- !
-                                          ! ---------------------------- !
+         ll_u_iterate = .TRUE.
+         ll_v_iterate = .TRUE.
+
+         DO i_inn = 1, nn_ninn_vp ! inner loop iterations
+         
+            ! mitgcm computes initial value of residual
+            jter = jter + 1
+            l_full_nf_update = jter == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
+
+            IF ( ll_u_iterate .OR. ll_v_iterate ) THEN
+            
+
+               zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate
+               zv_b(:,:) = v_ice(:,:)
+
+                                           ! ---------------------------- !
+               IF ( ll_u_iterate ) THEN    ! --- Solve for u-velocity --- !
+                                           ! ---------------------------- !
       
                   ! what follows could be subroutinized...
                   
-                  DO k = 1, nn_zebra_step ! "zebra" loop
-      
-                     IF ( k == 1 ); THEN; j_min = 2; ELSE; j_min = 3
-         
-                     zFU(:,:) = 0._wp
-         
-                     DO jj = j_min, jpjm1, nn_zebra_step
+                  DO jn = 1, nn_nzebra_vp ! "zebra" loop (! red-black-sor!!! )
+                  
+                     ! OPT: could be even better optimized with a true red-black SOR
+      
+                     IF ( jn == 1 ) THEN   ;   jj_min = 2 
+                     ELSE                  ;   jj_min = 3
+                     ENDIF
+         
+                     zFU(:,:)       = 0._wp
+                     zFU_prime(:,:) = 0._wp
+                     zBU_prime(:,:) = 0._wp
+         
+                     DO jj = jj_min, jpj - 1, nn_nzebra_vp
 
                         !-------------------------------------------
-                        ! 1) Tridiagonal system solver for each row
+                        ! -- Tridiagonal system solver for each row
                         !-------------------------------------------
+                        !
+                        ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan
+                        !
                         ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F
 
+                        ! - Right-hand side of tridiagonal system (zFU)
+                        DO ji = 2, jpi - 1    
+
+                           ! boundary condition substitution
+                           ! see Zhang and Hibler, 1997, Appendix B
+                           ! MV NOTE possibly not fully appropriate
+                           zAA3 = 0._wp
+                           IF ( ji == 2 )          zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj)
+                           IF ( ji == jpi - 1 )    zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj)
+     
+                           ! right hand side
+                           zFU(ji,jj) = ( zrhsu(ji,jj) &                                 ! right-hand side terms
+                               &      +   zAA3                                           ! boundary condition translation
+                               &      +   zDU(ji,jj) * u_ice(ji,jj-1)                    ! internal force, j-1
+                               &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1
+                    
+                        END DO
+                        
+                        ! - Thomas Algorithm 
+                        ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure)
+                        ! Forward sweep
+                        DO ji = 3, jpi - 1
+                           zw         = zAU(ji,jj) / zBU(ji-1,jj)
+                           zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj)
+                           zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj)
+                        END DO
+            
+                        ! Backward sweep
+                        ! MV I don't see how this will be reproducible
+                        u_ice(jpi-1,jj)     = zFU_prime(jpi-1,jj) / zBU_prime(jpi-1,jj) * umask(jpi-1,jj,1)                     ! do last row first
+                        DO ji = jpi-2 , 2, -1 ! all other rows    !   MV OPT: could be turned into forward loop (by substituting ji)
+                           u_ice(ji,jj)    = zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji,jj+1) * umask(ji,jj,1) / zBU_prime(ji,jj) !
+                        END DO            
+            
+                        !---------------       
+                        ! -- Relaxation
+                        !---------------
+                        DO ji = 2, jpi - 1    
+                           u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) )
+                        END DO
+
+                     END DO ! jj
+
+                  END DO ! zebra loop
+                  
+               ENDIF !   ll_u_iterate
+               
+               !                           ! ---------------------------- !
+               IF ( ll_v_iterate ) THEN    ! --- Solve for V-velocity --- !
+            !                              ! ---------------------------- !
+                                          
+                  ! MV OPT: what follows could be subroutinized...
+                  
+                  DO jn = 1, nn_nzebra_vp ! "zebra" loop
+      
+                     IF ( jn == 1 ) THEN   ;   ji_min = 2 
+                     ELSE                  ;   ji_min = 3
+                     ENDIF
+         
+                     zFV(:,:)       = 0._wp
+                     zFV_prime(:,:) = 0._wp
+                     zBV_prime(:,:) = 0._wp
+
+                     DO ji = ji_min, jpi - 1, nn_nzebra_vp 
+                        
+                        !!! It is intentional to have a ji then jj loop for V-velocity
+                        !!! ZH97 explain it is critical for convergence speed
+
+                        !-------------------------------------------
+                        ! -- Tridiagonal system solver for each row
+                        !-------------------------------------------
+                        ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F
+
                         ! --- Right-hand side of tridiagonal system (zFU)
-                        DO ji = fs_2, fs_jpim1    
-
-                           ! boundary condition substitution (check whether is correctly applied !!!)
+                        DO jj = 2, jpj - 1
+
+                           ! boundary condition substitution (check it is correctly applied !!!)
                            ! see Zhang and Hibler, 1997, Appendix B
                            zAA3 = 0._wp
-                           IF ( ji .EQ. fs_2 )     zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj)
-                           IF ( ji .EQ. fs_jpim1 ) zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj)
+                           IF ( jj .EQ. 2 )       zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1)
+                           IF ( jj .EQ. jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1)
      
                            ! right hand side
-                           zFU(ji,jj) = zrhsu(ji,jj) &               ! right-hand side terms
-                               &      + zAA3                        ! boundary condition translation
-                               &      + zDU(ji,jj) * u_ice(ji,jj-1) ! internal force, j-1
-                               &      + zEU(ji,jj) * u_ice(ji,jj+1) ! internal force, j+1
+                           zFV(ji,jj) = ( zrhsv(ji,jj) &              ! right-hand side terms
+                               &        + zAA3                        ! boundary condition translation
+                               &        + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1
+                               &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1
                     
-                           zFU(ji,jj) = zFU(ji,jj) * umask(ji,jj)     ! mask - useful ?
-        
                         END DO
        
@@ -777,136 +854,66 @@
                         ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure)
                         ! Forward sweep
-                        DO ji = fs_2 + 1, fs_jpim1
-                           zw         = zAU(ji,jj) / zBU(ji-1,jj)
-                           zBU(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj)
-                           zFU(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj)
+                        DO jj = 3, jpj - 1
+                           zw         = zAV(ji,jj) / zBV(ji,jj-1)
+                           zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1)
+                           zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 
                         END DO
             
                         ! Backward sweep
-                        u_ice(fs_jpim1,jj) = zFU(ji) / zBU(ji) ! last row
-                        DO ji  = fs_jpim - 1, fs_2, -1 ! can be turned into forward row by substituting ji if needed
-                           u_ice(ji,jj)    = zFU(ji) - zCU(ji,jj) * u_ice(ji,jj+1) / zBU(ji,jj)
+                        v_ice(ji,jpj-1)  = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1)  ! last row
+                        DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed
+                           v_ice(ji,jj)   = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj) / zBV_prime(ji,jj)
                         END DO            
             
                         !---------------       
-                        ! 2) Relaxation
+                        ! -- Relaxation
                         !---------------
-                        DO ji = fs_2, fs_jpim1    
-                           u_ice(ji,jj) = zu_p(ji,jj) + rn_relaxfacU * ( u_ice(ji,jj) - zu_p(ji,jj) )
+                        DO jj = 2, jpj - 1
+                           v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) )
                         END DO
 
-                     END DO ! jj
+                     END DO ! ji
 
                   END DO ! zebra loop
-                  
-                  CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1. )
-                  
-               ENDIF ! doIterate4u
-
-                                          ! ---------------------------- !
-               IF ( doIterate4v ) THEN    ! --- Solve for V-velocity --- !
-                                          ! ---------------------------- !
-                                          
-                  ! what follows could be subroutinized...
-                  
-                  DO k = 1, nn_zebra_step ! "zebra" loop
-      
-                     IF ( k == 1 ); THEN; i_min = 2; ELSE; i_min = 3
-         
-                     zFV(:,:) = 0._wp
-
-                     DO ji = i_min, fs_jpim1, nn_zebra_step 
-                        
-                        !!! HERE WE HAVE ji then jj loops - I don't know what this implies
-                        !!! But ZH explain it is critical for convergence speed
-
-                        !-------------------------------------------
-                        ! 1) Tridiagonal system solver for each row
-                        !-------------------------------------------
-                        ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F
-
-                        ! --- Right-hand side of tridiagonal system (zFU)
-                        DO jj = 2, jpjm1
-
-                           ! boundary condition substitution (check it is correctly applied !!!)
-                           ! see Zhang and Hibler, 1997, Appendix B
-                           zAA3 = 0._wp
-                           IF ( jj .EQ. 2 )     zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1)
-                           IF ( jj .EQ. jpjm1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1)
-     
-                           ! right hand side
-                           zFV(ji,jj) = zrhsv(ji,jj) &              ! right-hand side terms
-                               &      + zAA3                        ! boundary condition translation
-                               &      + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1
-                               &      + zEV(ji,jj) * v_ice(ji+1,jj) ! internal force, j+1
-                    
-                           zFV(ji,jj) = zFV(ji,jj) * vmask(ji,jj)     ! mask - useful ?
-        
-                        END DO
-       
-                        ! --- Thomas Algorithm 
-                        ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure)
-                        ! Forward sweep
-                        DO jj = 3, jpjm1
-                           zw         = zAV(ji,jj) / zBV(ji,jj-1)
-                           zBV(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1)
-                           zFV(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 
-                        END DO
-            
-                        ! Backward sweep
-                        v_ice(ji,jpjm1)  = zFV(ji) / zBV(ji) ! last row
-                        DO jj  = jpjm1 - 1, 2, -1 ! can be turned into forward row by substituting jj if needed
-                           v_ice(ji,jj)  = zFV(ji) - zCV(ji,jj) * v_ice(ji,jj+1) / zBV(ji,jj)
-                        END DO            
-            
-                        !---------------       
-                        ! 2) Relaxation
-                        !---------------
-                        DO jj = 2, jpjm1
-                           v_ice(ji,jj) = zv_p(ji,jj) + rn_relaxfacV * ( v_ice(ji,jj) - zv_p(ji,jj) )
-                        END DO
-
-                     END DO ! ji
-
-                  END DO ! zebra loop
-                  
-                  CALL lbc_lnk( 'icedyn_rhg_vp', v_ice, 'V', -1. )                                          
-                  
-               ENDIF ! doIterate4v
+                                    
+               ENDIF !   ll_v_iterate
+               
+               IF ( ( ll_u_iterate .OR. ll_v_iterate ) .OR. jter == nn_nvp ) CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
                               
-               !--------------------------------------------------
-               ! Continue the loop - based on velocity difference
-               !--------------------------------------------------
+               !--------------------------------------------------------------------------------------
+               ! -- Check convergence based on maximum velocity difference, continue or stop the loop
+               !--------------------------------------------------------------------------------------
 
                !------
                ! on U
                !------
-               IF ( doIterate4u .AND. MOD ( i_inner, nn_lsr_cvgcheck ) == 0 ) THEN
-               
-                  zuerr(:,:)   = 0._wp   
-                  DO jj = 2, jpjm1
-                     DO ji = 2, jpim1
-                        zuerr(ji,jj) = ( u_ice(ji,jj) - zu_p(ji,jj) ) * umask(ji,jj)
+               ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check
+               ! then we must optimize the use of the mpp_max, which is prohibitive                            
+                               
+               IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
+
+                  ! - Maximum U-velocity difference               
+                  zuerr(:,:) = 0._wp
+                  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
                      END DO
                   END DO
-                  
                   zuerr_max = MAXVAL( zuerr )
-                  CALL mpp_max( 'icedyn_rhg_evp', zuerr_max )   ! max over the global domain
-                  
-                  
-                  IF ( ln_rhg_chkcvg ) zdiag_u_idif(i_iter) = zuerr_max
-                   
-                  ! "safeguard against bad forcing" -> stop relaxation                    
-                  IF ( ( i_inner > 1 ) && zuerr_max > rn_uerr_max ) THEN
-                      !---> multiply relaxation factor above to stop relaxation
-                      zrelaxation_switchU = 0._wp
+                  CALL mpp_max( 'icedyn_rhg_evp', zuerr_max )   ! max over the global domain - damned!
+                  
+                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
+                  IF ( i_inn > 1 && zuerr_max > rn_uerr_max_vp ) THEN
+                      IF ( lwp ) " VP rheology error was too large : ", zu_err_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
+                      ll_u_iterate = .FALSE.
                   ENDIF
                   
-                  IF ( zuerr_max .LT. rn_uerr_min ) THEN                                        
-                      doIterate4u = .FALSE.
-                      IF ( lwp ) " VP rheology finished in outer U-iteration ", i_outer, " after ", i_inner, " iterations "
+                  ! - Stop if error small enough
+                  IF ( zuerr_max < rn_uerr_min_vp ) THEN                                        
+                      IF ( lwp ) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! "
+                      ll_u_iterate = .FALSE.
                   ENDIF
-                  
-               ENDIF
+                                               
+               ENDIF ! ll_u_iterate
 
                !------
@@ -914,69 +921,62 @@
                !------
                
-               IF ( doIterate4v .AND. MOD ( i_inner, nn_lsr_cvgcheck ) == 0 ) THEN
-               
+               IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
+               
+                  ! - Maximum V-velocity difference
                   zverr(:,:)   = 0._wp   
-                  DO jj = 2, jpjm1
-                     DO ji = 2, jpim1
-                        zverr(ji,jj) = ( v_ice(ji,jj) - zv_p(ji,jj) ) * vmask(ji,jj)
+                  DO jj = 2, jpj - 1
+                     DO ji = 2, jpi - 1
+                        zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1)
                      END DO
                   END DO
                   
                   zverr_max = MAXVAL( zverr )
-                  CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain
-                                   
-                  IF ( ln_rhg_chkcvg ) zdiag_v_idif(i_iter) = zverr_max
-                   
-                  ! "safeguard against bad forcing" -> stop relaxation                    
-                  IF ( ( i_inner > 1 ) && zverr_max > rn_uerr_max ) THEN
-                      !---> multiply relaxation factor above to stop relaxation
-                      zrelaxation_switchV = 0._wp
+                  CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain - damned!
+                  
+                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
+                  IF ( i_inn > 1 && zverr_max > rn_uerr_max_vp ) THEN
+                      IF ( lwp ) " VP rheology error was too large : ", zv_err_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
+                      ll_v_iterate = .FALSE.
                   ENDIF
                   
-                  IF ( zverr_max .LT. rn_verr_min ) THEN                                        
-                      doIterate4v = .FALSE.
-                      IF ( lwp ) " VP rheology finished in outer V-iteration ", i_outer, " after ", i_inner, " iterations "
+                  ! - Stop if error small enough
+                  IF ( zverr_max < rn_verr_min ) THEN                                        
+                      IF ( lwp ) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! "
+                      ll_v_iterate = .FALSE.
                   ENDIF
                   
-               ENDIF
-               
-               !-----------------------------------------------------------------------------
-               ! Convergence diagnostics: residual norm and mean kinetic energy at inner iteration
-               !-----------------------------------------------------------------------------
-               
-               IF( ln_rhg_chkcvg ) THEN
-               
-                  ! residual norm, velocity difference 
-                  ! ... zdiag_resnor
-                         
-                  ! mean kinetic energy            
-                  ztotke = 0._wp
-                  DO jj = 2, jpjm1
-                     DO ji = 2, jpim1
-                        ztotke = ztotke + ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * at_i(ji,jj) &
-                           &              ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * z1_32    
-                     END DO
-                  END DO
-            
-                  zdiag_meanke(i_iter) = ztotke / zglob_area
-                  
-               ENDIF
-            
-            ENDIF ! ---    end doIterate4u or doIterate4v
+               ENDIF ! ll_v_iterate
+               
+               !---------------------------------------------------------------------------------------
+               !
+               ! --- Calculate extra convergence diagnostics and save them
+               !
+               !---------------------------------------------------------------------------------------
+               IF( ln_rhg_chkcvg ) 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 )
+               
+            
+            ENDIF ! ---    end ll_u_iterate or ll_v_iterate
             
          END DO ! i_inn, end of inner loop
 
       !------------------------------------------------------------------------------!
+      !
       ! 6) Mask final velocities
+      !
       !------------------------------------------------------------------------------!
 
       END ! End of outer loop (i_out) =============================================================================================
 
-      !
       !------------------------------------------------------------------------------!
-      ! 7) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !
+      ! --- Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !
       !------------------------------------------------------------------------------!
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
+      !
+      ! OPT: subroutinize ?
+            
+      DO jj = 1, jpj - 1
+         DO ji = 1, jpi - 1
 
             ! shear at F points
@@ -988,6 +988,6 @@
       END DO           
       
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1 ! no vector loop
+      DO jj = 2, jpj - 1
+         DO ji = 2, jpi - 1 ! 
             
             ! tension**2 at T points
@@ -1017,8 +1017,12 @@
          END DO
       END DO
+      
       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 --- !
-      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+      ! 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 (:,:)
@@ -1027,7 +1031,11 @@
 
       !------------------------------------------------------------------------------!
-      ! 8) diagnostics
+      !
+      ! --- Diagnostics
+      !
       !------------------------------------------------------------------------------!
-      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
+      ! MV OPT: subroutinize
+      !
+      ! --- 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. &
          & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
@@ -1046,16 +1054,17 @@
       ENDIF
        
-      ! --- divergence, shear and strength --- !
+      ! --- Divergence, shear and strength --- !
       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('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
 
-      ! --- stress tensor --- !
+      ! --- 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) )
          !         
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
+         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) )  &
@@ -1066,13 +1075,11 @@
                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_vp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
          !
@@ -1086,10 +1093,12 @@
 
          DEALLOCATE( zsig1 , zsig2 , zsig3 )
+         
       ENDIF
 
-      ! --- SIMIP --- !
+      ! --- 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 !!!!!!!!!!!!!!!!
          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. )
@@ -1109,6 +1118,6 @@
             &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
          !
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
                ! 2D ice mass, snow mass, area transport arrays (X, Y)
                zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
@@ -1143,48 +1152,154 @@
       ENDIF
 
-      !
-      ! --- convergence diagnostics --- !
+      ! --- Convergence diagnostics --- !
       IF( ln_rhg_chkcvg ) THEN
           
-          ! still miss residual norm
-               
-         IF ( iom_use('mean_tke_ice') ) THEN
-            CALL iom_put( 'mean_tke_ice', zdiag_meanke )                  
-            DEALLOCATE( zdiag_meanke )       
-         ENDIF
-         
-         IF ( iom_use('u_veldif') ) THEN
-            CALL iom_put ( zdiag_u_idif(i_iter) )
-            DEALLOCATE( zdiag_u_idif )
-         ENDIF
-         
-         IF ( iom_use('v_veldif') ) THEN
-            CALL iom_put ( zdiag_v_idif(i_iter) )
-            DEALLOCATE( zdiag_v_idif )
-         ENDIF
-         
+         IF( iom_use('uice_cvg') ) THEN
+            CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , &
+                  &                        ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
+         ENDIF   
+        
       ENDIF ! ln_rhg_chkcvg
 
-! EVP convergence diag      
-!         IF( iom_use('uice_cvg') ) THEN
-!         
-!        
-!           ! EVP-diag
-!           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_vp
    
+   
+   
+   SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, &
+                  &       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  ***
+      !!                     
+      !! ** Purpose :   check convergence of VP ice rheology
+      !!
+      !! ** Method  :   create a file ice_cvg.nc containing a few convergence diagnostics
+      !!                This routine is called every sub-iteration, so it is cpu expensive
+      !!
+      !!                Calculates / stores
+      !!                   - maximum absolute U-V difference (uice_cvg, u_dif, v_dif, m/s)
+      !!                   - residuals in U, V and UV-mean taken as square-root of area-weighted mean square residual (u_res, v_res, vel_res, N/m2)
+      !!                   - mean kinetic energy (mke_ice, J/m2)
+      !!                   
+      !!
+      !! ** 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, pmt              ! now velocity and mass per unit area 
+      REAL(wp),                 INTENT(in) ::   puerr_max, pverr_max     ! absolute mean velocity difference
+      REAL(wp),                 INTENT(in) ::   pglob_area               ! global ice area
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients 
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV
+      !!
+      INTEGER           ::   it, idtime, istatus
+      INTEGER           ::   ji, jj          ! dummy loop indices
+      REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zmke, zu, zv ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj)   ::   zu_res(:,:), zv_res(:,:), zvel2(:,:) ! local arrays
+                                                                             
+      CHARACTER(len=20) ::   clname
+      ::   zres           ! check convergence
+      !!----------------------------------------------------------------------
+
+      ! create file
+      IF( kt == nit000 .AND. kiter == 1 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'rhg_cvg_vp : 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_ucvg )
+            istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE   , (/ idtime /), nvarid_ures )
+            istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE   , (/ idtime /), nvarid_vres )
+            istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE   , (/ idtime /), nvarid_velres )
+            istatus = NF90_DEF_VAR( ncvgid, 'u_dif', NF90_DOUBLE   , (/ idtime /), nvarid_udif )
+            istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE   , (/ idtime /), nvarid_vdif )
+            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE   , (/ idtime /), nvarid_mke )
+            istatus = NF90_ENDDEF(ncvgid)
+         ENDIF
+         !
+      ENDIF
+
+      ! time
+      it = ( kt - 1 ) * kitermax + kiter
+
+      ! --- Max absolute velocity difference with previous iterate (zveldif)
+      ! EVP code      
+!      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
+!         zveldif = 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
+!         zveldif = MAXVAL( zres )
+!         CALL mpp_max( 'icedyn_rhg_vp', zveldif )   ! max over the global domain
+!      ENDIF
+      ! VP code
+      zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code 
+                                            ! if puerrmask and pverrmax are masked at 15% (TEST)
+      
+      ! -- Mean residual (N/m^2), zu_res_mean
+      ! Here we take the residual of the linear system (N/m^2), 
+      ! We define it as in mitgcm: square-root of area-weighted mean square residual
+      ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 
+      ! i.e., how close we are to a solution
+      DO jj = 2, jpj - 1
+         DO ji = 2, jpi - 1                                      
+            zu_res(ji,jj)  = zrhsu(ji,jj) + zDU(ji,jj) * pu(ji,jj-1) + zEU(ji,jj) * pu(ji,jj+1) &
+               &           - zAU(ji,jj) * pu(ji-1,jj) - zBU(ji,jj) * pu(ji,jj) - zCU(ji,jj) * pu(ji+1,jj)
+                           
+            zv_res(ji,jj)  = zrhsv(ji,jj) + zDV(ji,jj) * pv(ji-1,jj) + zEV(ji,jj) * pv(ji+1,jj) &
+               &           - zAV(ji,jj) * pv(ji,jj-1) - zBV(ji,jj) * pv(ji,jj) - zCV(ji,jj) * pv(ji,jj+1)                           
+          END DO
+       END DO                  
+       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:) )
+       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:) )
+       zu_res_mean = SQRT( zu_resmean / pglob_area )
+       zv_res_mean = SQRT( zv_resmean / pglob_area )
+       zvelres     = MEAN( zu_res_mean, zv_res_mean )
+                                          
+       ! -- Global mean kinetic energy per unit area (J/m2)  
+       DO jj = 2, jpj - 1
+          DO ji = 2, jpi - 1                   
+             zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point
+             zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) )
+             zvel2(:,:)  = zu*zu + zv*zv              ! square of ice velocity at T-point  
+          END DO
+       END DO
+       
+       zmke = 0.5_wp * globsum( 'ice_rhg_vp', zmt(:,:) * e1e2t(:,:) * zvel2(:,:) )
+                  
+         !                                                ! ==================== !
+
+      IF( lwm ) THEN
+         ! write variables
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) )
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) )
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) )
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) )
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) )
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) )
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) )
+         ! close file
+         IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid )
+      ENDIF
+      
+   END SUBROUTINE rhg_cvg_vp
+   
+
 
    SUBROUTINE rhg_vp_rst( cdrw, kt )
