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 13450)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90	(revision 13493)
@@ -34,7 +34,9 @@
    INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology
 !! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
+   INTEGER, PARAMETER ::   np_rhgVP  = 3   ! VP rheology
 
    ! ** namelist (namrhg) **
    LOGICAL ::   ln_rhg_EVP       ! EVP rheology
+   LOGICAL ::   ln_rhg_VP        ! EVP rheology
    !
    !! * Substitutions
@@ -76,8 +78,15 @@
       !--------------!   
       SELECT CASE( nice_rhg )
-      !                                !------------------------!
-      CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic !
-         !                             !------------------------!
+      !                                !---------------------------------!
+      CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic rheology !
+         !                             !---------------------------------!
          CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
+         !         
+      END SELECT
+
+      !                                !---------------------------------------------!
+      CASE( np_rhgVP  )                ! Viscous-Plastic rheology                    !
+         !                             !---------------------------------------------!
+         CALL ice_dyn_rhg_vp ( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
          !         
       END SELECT
@@ -85,4 +94,5 @@
       IF( lrst_ice ) THEN                       !* write EVP fields in the restart file
          IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt )
+         IF( ln_rhg_VP  )   CALL rhg_lsr_rst( 'WRITE', kt )
       ENDIF
       !
@@ -110,5 +120,6 @@
       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
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg, &
+    &                       ln_rhg_VP 
       !!-------------------------------------------------------------------
       !
@@ -133,4 +144,6 @@
          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,*) ''    
+         WRITE(numout,*) '      rheology VP     (icedyn_rhg_lsr)                     ln_rhg_VP     = ', ln_rhg_VP
       ENDIF
       !
@@ -139,7 +152,9 @@
       IF( ln_rhg_EVP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEVP    ;   ENDIF
 !!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF
+      IF( ln_rhg_VP  ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgVP     ;   ENDIF
       IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' )
       !
       IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
+      IF( ln_rhg_VP   )   CALL rhg_lsr_rst( 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE ice_dyn_rhg_init
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 13493)
@@ -0,0 +1,1254 @@
+! TODOLIST
+
+!
+! Code residual norm of the linear system
+! Verify convergence check
+!
+! Define all symbols
+! Check uses
+! 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
+!
+! Subroutinize tridiagonal systems ?
+! 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
+
+MODULE icedyn_rhg_vp
+   !!======================================================================
+   !!                     ***  MODULE  icedyn_rhg_vp  ***
+   !!   Sea-Ice dynamics : Viscous-plastic rheology with LSR technique
+   !!======================================================================
+   !!
+   !! History :   -   !  1997-20  (J. Zhang, M. Losch) Original code, implementation into mitGCM
+   !!            4.0  !  2020-09  (M. Vancoppenolle) Adaptation to SI3
+   !!
+   !!----------------------------------------------------------------------
+#if defined key_si3
+   !!----------------------------------------------------------------------
+   !!   'key_si3'                                       SI3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_dyn_rhg_vp : computes ice velocities from VP rheolog with LSR solvery
+   !!   rhg_vp_rst     : read/write EVP fields in ice restart
+   !!----------------------------------------------------------------------
+   USE phycst         ! Physical constants
+   USE dom_oce        ! Ocean domain
+   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
+   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
+   USE ice            ! sea-ice: ice variables
+   USE icevar         ! ice_var_sshdyn
+   USE icedyn_rdgrft  ! sea-ice: ice strength
+   USE bdy_oce , ONLY : ln_bdy 
+   USE bdyice 
+#if defined key_agrif
+   USE agrif_ice_interp
+#endif
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE lbclnk         ! lateral boundary conditions (or mpp links)
+   USE prtctl         ! Print control
+
+   USE netcdf         ! NetCDF library for convergence test
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_dyn_rhg_vp   ! called by icedyn_rhg.F90
+   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
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
+   !! $Id: icedyn_rhg_vp.F90 13279 2020-07-09 10:39:43Z clem $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_dyn_rhg_vp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
+      !!-------------------------------------------------------------------
+      !!
+      !!                 ***  SUBROUTINE ice_dyn_rhg_vp  ***
+      !!                             VP-LSR-C-grid
+      !!
+      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
+      !!  stress and sea-surface slope. Internal forces assume viscous-plastic rheology (Hibler, 1979)
+      !!  
+      !!  The resolution algorithm  follows from Zhang and Hibler (1997) and Losch (2010)
+      !!  
+      !!  It is based on a linear slope relaxation algorithm (LSR)
+      !!
+      !!  * Time stepping
+      !!
+      !!  The idea is to arrange the momentum equations as follows:
+      !!  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 )
+      !!  
+      !!  * Spatial discretization 
+      !!
+      !!  Assumes a C-grid
+      !!
+      !!  The points in the C-grid look like this, my darling
+      !!
+      !!                              (ji,jj)
+      !!                                 |
+      !!                                 |
+      !!                      (ji-1,jj)  |  (ji,jj)
+      !!                             ---------   
+      !!                            |         |
+      !!                            | (ji,jj) |------(ji,jj)
+      !!                            |         |
+      !!                             ---------   
+      !!                     (ji-1,jj-1)     (ji,jj-1)
+      !!
+      !! ** Inputs  : - wind forcing (stress), oceanic currents
+      !!                ice total volume (vt_i) per unit area
+      !!                snow total volume (vt_s) per unit area
+      !!
+      !! ** Action  : 
+      !!             
+      !! ** Steps   :
+      !!            
+      !! ** Notes   : 
+      !!             
+      !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010. 
+      !!             
+      !!             
+      !!-------------------------------------------------------------------
+      !!
+!!!      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
+!!!      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), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)      
+                  
+      !! --- 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
+      !------------------------------------------------------------------------------!
+      zrhoco = rau0 * rn_cio 
+
+      ! ecc2: square of yield ellipse eccentricity
+      ecc2    = rn_ecc * rn_ecc
+      z1_ecc2 = 1._wp / ecc2
+         
+      ! Initialise stress tensor from previous time step
+      zs1 (:,:) = pstress1_i (:,:) 
+      zs2 (:,:) = pstress2_i (:,:)
+      zs12(:,:) = pstress12_i(:,:)
+      
+      ! Convergence checks
+      IF( ln_rhg_chkcvg ) THEN 
+      
+         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)
+      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
+      ELSE                         ;   zkt = 0._wp
+      ENDIF
+
+      !------------------------------------------------------------------------------!
+      ! 2) Time-independent quantities
+      !------------------------------------------------------------------------------!
+      
+      CALL ice_strength ! strength at T points
+      
+      !--------
+      ! F-mask      
+      !--------
+      ! 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.
+            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
+      END DO
+      CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp )
+
+      ! Lateral boundary conditions on velocity (modify zfmask)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            IF( zfmask(ji,jj) == 0._wp ) THEN
+               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
+                  &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
+            ENDIF
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         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) ) )
+        ENDIF
+      END DO
+      DO ji = 2, jpim1
+         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) ) )
+         ENDIF
+      END DO
+      
+      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
+
+            ! Ice fraction at U-V points
+            zaU(ji,jj)      = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zaV(ji,jj)      = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+
+            ! Snow and ice mass at U-V points
+            zm1             = ( 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)
+                          
+            ! Mass per unit time
+            zmassU_t(ji,jj) = zmassU * r1_rdtice
+            zmassV_t(ji,jj) = zmassV * r1_rdtice
+
+            ! Acceleration term contribution to RHS (depends on velocity at previous time step)            
+            zmU_t(ji,jj)    = zmassU_t(ji,jj) * u_ice(ji,jj)
+            zmV_t(ji,jj)    = zmassV_t(ji,jj) * v_ice(ji,jj)
+            
+            ! Ocean currents at U-V points
+            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
+            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
+
+            ! Coriolis factor at T points (m*f)
+            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
+            
+            ! Wind stress
+            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+
+            ! Force due to sea surface tilt(- m*g*GRAD(ssh))
+            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! masks
+            zmsk00x(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
+            ! switches
+            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
+
+         END DO
+      END DO   
+            
+      !------------------------------------------------------------------------------!
+      ! 4) Start outer loop, update velocities
+      !------------------------------------------------------------------------------!
+      
+      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
+      !------------------------------------------------------------------------------!
+
+         !      
+         ! 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
+         !--------------------------------------------------------
+
+         ! --- 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
+
+               ! 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)   &
+                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+            END DO
+         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
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+              
+               ! 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)   &
+                  &    ) * 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)   &
+                  &   ) * r1_e1e2t(ji,jj)
+               zdt2 = zdt * zdt
+               
+               ! delta at T points
+               zdeltat = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+               
+               ! delta* at T points (following Lemieux and Dupont, GMD 2020)
+               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
+               
+               ! P/delta at T-points
+               zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltat_star
+               
+               ! Temporary zzt and zet factors at T-points
+               zzt(ji,jj)     = zp_deltastar_t * r1_e1e2t(ji,jj)
+               zet(ji,jj)     = zzt(ji,jj)     * z1_ecc2 
+                          
+            END DO
+         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
+               
+               ! P/delta* at F points
+               zp_deltastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) )
+               
+               ! Temporary zef factor at F-point
+               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2
+
+            END DO
+         END DO
+         
+         CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. )
+         
+         !---------------------------------------------------
+         ! 5.2 Ocean-ice drag and Coriolis RHS contributions
+         !---------------------------------------------------
+         
+         DO jj = 2, jpjm1
+             DO ji = fs_2, fs_jpim1
+                
+                !--- 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)
+
+                !--- 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) ) )
+                 
+                !--- Ocean-ice drag contributions to RHS 
+                ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj)
+                ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj)
+                
+                ! --- U-component of Coriolis Force (energy conserving formulation)
+                ! 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) ) )
+                           
+                ! --- 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) ) )
+         
+             END DO
+         END DO
+         
+         ! ??? lbclnk ???
+
+         !--------------------------------------
+         ! 5.3 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
+            
+               ! 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 )
+                                            
+               ! 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) 
+
+               ! 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 )
+
+               ! 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)
+
+            END DO
+         END DO
+         
+         ! ??? lbclnk ???
+         
+         ! --- 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
+         ! My guess is that this is the way to enforce boundary conditions on strain rate tensor
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               
+               ! 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)
+               
+               ! 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)
+
+            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               
+               ! --- U component of internal force contribution to RHS at U points
+               zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & 
+                              (    e2u(ji,jj)    * ( zs1_rhsu(ji+1,jj) - zs1_rhsu(ji,jj) )                                                                 &
+                  &      +         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) * &
+                  &           (    e1v(ji,jj)    * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) )                                                                 &
+                  &      +         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
+         !---------------------------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+            
+               ! still miss ice ocean stress and acceleration contribution
+               zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj)
+               zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsu(ji,jj)
+
+            END DO
+         END DO
+         
+      ! ??? lbclnk ???
+     
+      !---------------------------------------------------------------------------------------!
+      ! 6) Linear system matrix
+      !---------------------------------------------------------------------------------------!
+      
+         ! Linear system matrix contains all implicit contributions 
+         ! 1) internal forces, 2) acceleration, 3) ice-ocean drag
+
+         ! The linear system equation is written as follows
+         ! 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 
+         !         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) 
+               
+      !------------------------------------------------------------------------------!
+      ! 6) 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 --- !
+                                          ! ---------------------------- !
+      
+                  ! 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
+
+                        !-------------------------------------------
+                        ! 1) Tridiagonal system solver for each row
+                        !-------------------------------------------
+                        ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F
+
+                        ! --- Right-hand side of tridiagonal system (zFU)
+                        DO ji = fs_2, fs_jpim1    
+
+                           ! boundary condition substitution (check whether 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)
+     
+                           ! 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
+                    
+                           zFU(ji,jj) = zFU(ji,jj) * umask(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 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)
+                        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)
+                        END DO            
+            
+                        !---------------       
+                        ! 2) 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) )
+                        END DO
+
+                     END DO ! jj
+
+                  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
+                              
+               !--------------------------------------------------
+               ! Continue the loop - based on velocity difference
+               !--------------------------------------------------
+
+               !------
+               ! 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)
+                     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
+                  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 "
+                  ENDIF
+                  
+               ENDIF
+
+               !------
+               ! on V
+               !------
+               
+               IF ( doIterate4v .AND. MOD ( i_inner, nn_lsr_cvgcheck ) == 0 ) THEN
+               
+                  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)
+                     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
+                  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 "
+                  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
+            
+         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) 
+      !------------------------------------------------------------------------------!
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+
+            ! shear at F points
+            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+         END DO
+      END DO           
+      
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! tension**2 at T points
+            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               &   ) * r1_e1e2t(ji,jj)
+            zdt2 = zdt * zdt
+            
+            ! shear**2 at T points (doc eq. A16)
+            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+            
+            ! shear at T points
+            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! divergence at T points
+            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+               &             ) * r1_e1e2t(ji,jj)
+            
+            ! delta at T points
+            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 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. )
+      pstress1_i (:,:) = zs1 (:,:)
+      pstress2_i (:,:) = zs2 (:,:)
+      pstress12_i(:,:) = zs12(:,:)
+      !
+
+      !------------------------------------------------------------------------------!
+      ! 8) diagnostics
+      !------------------------------------------------------------------------------!
+      ! --- 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
+
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
+            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+         !
+         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
+
+      ENDIF
+       
+      ! --- 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 --- !
+      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
+               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
+                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
+
+               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
+
+               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
+
+!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
+!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
+!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
+!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
+               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
+               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
+               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
+         !
+         CALL iom_put( 'isig1' , zsig1 )
+         CALL iom_put( 'isig2' , zsig2 )
+         CALL iom_put( 'isig3' , zsig3 )
+         !
+         ! Stress tensor invariants (normal and shear stress N/m)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
+
+         DEALLOCATE( zsig1 , zsig2 , zsig3 )
+      ENDIF
+
+      ! --- SIMIP --- !
+      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
+         !
+         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. )
+
+         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
+      ENDIF
+
+      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
+         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
+         !
+         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
+            &      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
+               ! 2D ice mass, snow mass, area transport arrays (X, Y)
+               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
+               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
+
+               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
+               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+
+               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+
+               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
+               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
+
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
+            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
+            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+
+         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
+         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport 
+         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
+         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
+         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
+         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
+
+         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
+            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
+
+      ENDIF
+
+      !
+      ! --- 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
+         
+      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_vp_rst( cdrw, kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rhg_vp_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write RHG file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+      !
+      INTEGER  ::   iter            ! local integer
+      INTEGER  ::   id1, id2, id3   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
+            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
+            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
+               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+            ELSE                                     ! start rheology from rest
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
+               stress1_i (:,:) = 0._wp
+               stress2_i (:,:) = 0._wp
+               stress12_i(:,:) = 0._wp
+            ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
+            stress1_i (:,:) = 0._wp
+            stress2_i (:,:) = 0._wp
+            stress12_i(:,:) = 0._wp
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
+         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+         !
+         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+         !
+      ENDIF
+      !
+   END SUBROUTINE rhg_vp_rst
+
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Empty module           NO SI3 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!==============================================================================
+END MODULE icedyn_rhg_vp
+
