Index: /NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90	(revision 13590)
+++ /NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90	(revision 13591)
@@ -138,4 +138,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
       !
@@ -170,5 +171,5 @@
       !! --- diags
       REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i         
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
       !! --- SIMIP diags
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
@@ -809,4 +810,5 @@
       IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
 
@@ -830,5 +832,5 @@
                ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
                zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure
-               zsig_II(ji,jj)   =   - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress
+               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress
                
             END DO
@@ -852,5 +854,5 @@
       IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
          !
-         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )         
+         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !         
          DO jj = 2, jpj - 1
@@ -867,10 +869,10 @@
                ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
                zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure
-               zsig_II(ji,jj)   =   - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress
+               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )     ! 2nd  ''       '', aka maximum shear stress
 
                ! Normalized  principal stresses (used to display the ellipse)
-               z1_strength      =   1._wp / strength(ji,jj)
-               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  - zsig_II(ji,jj) ) * z1_strength
-               zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+               z1_strength      =   1._wp / MAX ( 1.0, strength(ji,jj) )
+               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
             END DO
          END DO
@@ -882,5 +884,5 @@
          CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
 
-         DEALLOCATE( zsig1_p , zsig2_p )
+         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
          
       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 13590)
+++ /NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13591)
@@ -50,4 +50,6 @@
    INTEGER ::   nvarid_vdif
    INTEGER ::   nvarid_mke
+   INTEGER ::   nvarid_usub, nvarid_vsub
+
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
 
@@ -129,8 +131,11 @@
       LOGICAL ::   ll_u_iterate, ll_v_iterate   ! continue iteration or not
       !
-      INTEGER ::   ji, jj, jn          ! dummy loop indices
+      INTEGER ::   ji, ji2, 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_thomas           ! tridiagonal system version (1=mitgcm, 2=martin)
+
       INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)      
       !
@@ -146,5 +151,5 @@
       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) ::   zAA3, zw, ztau, zuerr_max, zverr_max
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
@@ -154,4 +159,5 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension
       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
@@ -182,4 +188,5 @@
       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) ::   zCU_prime, zCV_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)
@@ -193,5 +200,5 @@
       !! --- diags
       REAL(wp) ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y
-      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f, zten_i   ! stress tensor components
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f           ! stress tensor components
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi
@@ -201,6 +208,18 @@
                         
       !!----------------------------------------------------------------------------------------------------------------------
+      ! DEBUG put all forcing terms to zero
+         ! air-ice drag
+         utau_ice(:,:) = 0._wp
+         vtau_ice(:,:) = 0._wp
+         ! coriolis
+         ff_t(:,:) = 0._wp
+         ! ice-ocean drag
+         rn_cio = 0._wp
+         ! ssh 
+         ! done line 330 !!! dont forget to act there
+      ! END DEBUG
 
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)'
+      IF( lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)'
             
       !------------------------------------------------------------------------------!
@@ -219,7 +238,12 @@
       END DO
       
-      IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2; ELSE; nn_zebra_vp = 1; ENDIF 
-      
+      IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2
+                         ELSE; nn_zebra_vp = 1; ENDIF 
+
       nn_nvp = nn_nout_vp * nn_ninn_vp ! maximum number of iterations
+
+      IF( lwp )   WRITE(numout,*) ' ln_zebra_vp : ', ln_zebra_vp
+      IF( lwp )   WRITE(numout,*) ' nn_zebra_vp : ', nn_zebra_vp
+      IF( lwp )   WRITE(numout,*) ' nn_nvp      : ', nn_nvp
       
       zrhoco = rau0 * rn_cio 
@@ -242,4 +266,10 @@
       ELSE                         ;   zkt = 0._wp
       ENDIF
+
+      zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp
+      zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp;
+      zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp;
+      zrhsu(:,:) = 0._wp; zrhsv(:,:) = 0._wp
+      zf_rhsu(:,:) = 0._wp; zf_rhsv(:,:) = 0._wp
 
       !------------------------------------------------------------------------------!
@@ -304,4 +334,5 @@
       !    non-embedded sea ice: use ocean surface for slope calculation
       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
+      zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!!
       
       DO jj = 2, jpj - 1
@@ -349,10 +380,20 @@
             ! Mask for lots of ice (1) or little ice (0)
             IF( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
-            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+            ELSE                                                     ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
             IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
-            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
+            ELSE                                                     ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
 
          END DO
       END DO   
+
+      CALL iom_put( 'zmsk00x'    , zmsk00x  )   ! MV DEBUG
+      CALL iom_put( 'zmsk00y'    , zmsk00y  )   ! MV DEBUG
+      CALL iom_put( 'zmsk01x'    , zmsk01x  )   ! MV DEBUG
+      CALL iom_put( 'zmsk01y'    , zmsk01y  )   ! MV DEBUG
+      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG
+      CALL iom_put( 'ztaux_ai'   , ztaux_ai )   ! MV DEBUG
+      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG
+      CALL iom_put( 'zspgU'      , zspgU    )   ! MV DEBUG
+      CALL iom_put( 'zspgV'      , zspgV    )   ! MV DEBUG
             
       !------------------------------------------------------------------------------!
@@ -368,4 +409,6 @@
 
       DO i_out = 1, nn_nout_vp
+
+         IF( lwp )   WRITE(numout,*) ' outer loop  i_out : ', i_out
                
          ! Velocities used in the non linear terms are the average of the past two iterates
@@ -400,6 +443,9 @@
             END DO
          END DO
-         
+
          CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan
+         CALL iom_put( 'zds'        , zds      )   ! MV DEBUG
+
+         IF( lwp )   WRITE(numout,*) ' outer loop  1a i_out : ', i_out
 
          DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
@@ -441,4 +487,10 @@
          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. )
 
+         CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG
+         CALL iom_put( 'zet'        , zet      )   ! MV DEBUG
+         CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG
+
+         IF( lwp )   WRITE(numout,*) ' outer loop  1b i_out : ', i_out
+
          DO jj = 1, jpj - 1
             DO ji = 1, jpi - 1
@@ -454,9 +506,11 @@
          
          CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. )
+         CALL iom_put( 'zef'           , zef            ) ! MV DEBUG
+         IF( lwp )   WRITE(numout,*) ' outer loop  1c i_out : ', i_out
 
          !---------------------------------------------------
          ! -- Ocean-ice drag and Coriolis RHS contributions
          !---------------------------------------------------
-         
+
          DO jj = 2, jpj - 1
              DO ji = 2, jpi - 1
@@ -489,4 +543,15 @@
              END DO
          END DO
+
+         IF( lwp )   WRITE(numout,*) ' outer loop  1d i_out : ', i_out
+         CALL lbc_lnk( 'icedyn_rhg_vp', zCwU, 'U', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_vp', zCwV, 'V', 1. )
+
+         CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG
+         CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG
+         IF( lwp )   WRITE(numout,*) ' outer loop  1e i_out : ', i_out
+         CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG
+         CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG
+         IF( lwp )   WRITE(numout,*) ' outer loop  1f i_out : ', i_out
          
          ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system, 
@@ -515,4 +580,9 @@
             END DO
          END DO
+
+         CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG
+         CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG
+         CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG
+         CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG
          
          ! a priori, no lbclnk, because rhsu is only used in the inner domain
@@ -522,4 +592,6 @@
          ! My guess is that this is the way to enforce boundary conditions on strain rate tensor
 
+         IF( lwp )   WRITE(numout,*) ' outer loop 2 i_out : ', i_out
+
          DO jj = 1, jpj - 1
             DO ji = 1, jpi - 1
@@ -533,5 +605,11 @@
             END DO
          END DO
-         
+
+         CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. )
+
+         CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG
+         CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG
+
          ! a priori, no lbclnk, because rhsu are only used in the inner domain
 
@@ -555,4 +633,7 @@
             END DO
          END DO
+
+         CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG
+         CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG
          
          !---------------------------
@@ -571,5 +652,12 @@
          END DO
          
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', 1., zrhsv, 'V',  1. )
+
+         CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG
+         CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG
+         
          ! inner domain calculations -> no lbclnk
+
+         IF( lwp )   WRITE(numout,*) ' outer loop 4 i_out : ', i_out
      
          !---------------------------------------------------------------------------------------!
@@ -662,5 +750,22 @@
             END DO
          END DO
-               
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU  , 'U', 1., zAV  , 'V',  1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU  , 'U', 1., zBV  , 'V',  1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU  , 'U', 1., zCV  , 'V',  1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zDU  , 'U', 1., zDV  , 'V',  1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zEU  , 'U', 1., zEV  , 'V',  1. )
+               
+         CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG
+         CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG
+         CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG
+         CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG
+         CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG
+         CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG
+         CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG
+         CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG
+         CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG
+         CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG
+
       !------------------------------------------------------------------------------!
       !
@@ -673,15 +778,27 @@
          ll_v_iterate = .TRUE.
 
+         ! zBU
+         ! it is probably the way zBU is used that is problematic...
+         ! because zBU itself does not have the zipper!!!
+         ! BUG only occurs when zBU is different from zero
+         ! H1 - the jj loop does not work well with jpj = 149 ?
+         ! zBU(:,:) = 0._wp; zBV(:,:) = 0._wp
+
          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( lwp )   WRITE(numout,*) ' inner loop 1 i_inn : ', i_inn
+         
+            !--- mitgcm computes initial value of residual here...
+
+            jter             = jter + 1
+            ! l_full_nf_update = jter == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
+
+            zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate
+            zv_b(:,:) = v_ice(:,:)
+
+            zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp
+            zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp
+
+            IF ( ll_u_iterate .OR. ll_v_iterate )   THEN
 
                                            ! ---------------------------- !
@@ -689,4 +806,5 @@
                                            ! ---------------------------- !
       
+                  ! Thomas Algorithm  for tridiagonal solver
                   ! what follows could be subroutinized...
                   
@@ -698,14 +816,13 @@
                      ELSE                  ;   jj_min = 3
                      ENDIF
-         
-                     zFU(:,:)       = 0._wp
-                     zFU_prime(:,:) = 0._wp
-                     zBU_prime(:,:) = 0._wp
-         
+
+                     IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min
+
+                     !----------------------------------------------------------
+                     ! -- Boundary condition (link with neighbouring processor)
+                     !----------------------------------------------------------
+
                      DO jj = jj_min, jpj - 1, nn_zebra_vp
 
-                        !-------------------------------------------
-                        ! -- Tridiagonal system solver for each row
-                        !-------------------------------------------
                         !
                         ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan
@@ -722,32 +839,110 @@
                            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
-                    
+                           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
+
+                     END DO
+
+                     ! CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. )
                         
-                        ! - 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)
+                     !---------------
+                     ! Forward sweep
+                     !---------------
+                     nn_thomas = 2 ! 1 = martin, 2 = mitgcm
+
+                     IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm
+
+                        DO jj = jj_min, jpj - 1, nn_zebra_vp
+
+                           DO ji = 3, jpi - 1
+
+                              zfac             = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) )
+                              zw               = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 ) 
+                              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
+                        
                         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
-                        !---------------
+
+                     ELSE ! nn_thomas == 2      ! mitGCM algorithm
+
+                        DO jj = jj_min, jpj - 1, nn_zebra_vp
+
+                           ! ji = 2
+                           zw                   =   zBU(2,jj)
+                           zfac                 =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) )
+                           zw                   =   zfac * 1.0_wp /  MAX ( ABS( zw ) , epsi20 )
+                           zCU_prime(2,jj)      =   zw * zCU(2,jj)
+                           zFU_prime(2,jj)      =   zw * zFU(2,jj)
+
+                           DO ji = 3, jpi - 1
+
+                              zw                =   zBU(ji,jj) - zAU(ji,jj) * zCU_prime(ji-1,jj)
+                              zfac              =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) )
+                              zw                =   zfac * 1.0_wp /  MAX ( ABS( zw ) , epsi20 )
+                              zCU_prime(ji,jj)  =   zw * zCU(ji,jj) 
+                              zFU_prime(ji,jj)  =   zw * ( zFU(ji,jj) - zAU(ji,jj) * zFU_prime(ji-1,jj) )
+
+                           END DO
+
+                        END DO
+
+                     ENDIF
+
+                     ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zFU_prime, 'U', 1. )
+
+                     !----------------
+                     ! Backward sweep
+                     !----------------
+                     IF ( nn_thomas == 1 ) THEN ! MV version
+
+                        DO jj = jj_min, jpj - 1, nn_zebra_vp
+
+                           ! last row 
+                           zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) )
+                           u_ice(jpi-1,jj)     = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 
+                                  &            * umask(jpi-1,jj,1)
+
+                           DO ji2 = 2 , jpi-2 ! all other rows       
+                          !DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop
+                              ji =  jpi - ji2
+                              ! ji2 = 2       -> ji = jpi - 2
+                              ! ji2 = jpi - 2 -> ji = 2
+                              zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) )
+                              u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) * umask(ji,jj,1) ) & ! this
+                                   ! line is guilty
+                                              &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) 
+                           END DO            
+
+                        END DO
+
+                     ELSE ! nn_thomas == 2 ! mitGCM version
+
+                        DO jj = jj_min, jpj - 1, nn_zebra_vp
+
+                           u_ice(jpi-1,jpj)   =   zFU_prime(jpi-1,jpj)
+
+                           DO ji2 = 2 , jpi-2 ! all other rows       !   MV OPT: could be turned into forward loop (by substituting ji)
+                              ji =  jpi - ji2
+                              u_ice(ji,jj)    =   zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj)
+                           END DO
+                     
+                        END DO
+
+                     ENDIF
+
+                     !---------------       
+                     ! -- Relaxation
+                     !---------------
+
+                     DO jj = jj_min, jpj - 1, nn_zebra_vp
+            
                         DO ji = 2, jpi - 1    
                         
@@ -755,20 +950,19 @@
                            
                            ! velocity masking for little-ice and no-ice cases
-                           ! put 0.01 of ocean current, appropriate choice to avoid too much spreading of the ice edge
                            u_ice(ji,jj) =   zmsk00x(ji,jj)                                        & 
                               &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             &
-                                          + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    )
+                              &           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) * umask(ji,jj,1)
                            
                         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...
@@ -779,9 +973,7 @@
                      ELSE                  ;   ji_min = 3
                      ENDIF
-         
-                     zFV(:,:)       = 0._wp
-                     zFV_prime(:,:) = 0._wp
-                     zBV_prime(:,:) = 0._wp
-
+
+                     IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min
+         
                      DO ji = ji_min, jpi - 1, nn_zebra_vp 
                         
@@ -800,6 +992,6 @@
                            ! 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. jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1)
+                           IF ( jj == 2 )       zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1)
+                           IF ( jj == jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1)
      
                            ! right hand side
@@ -808,12 +1000,13 @@
                                &        + 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
-                    
+
                         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, jpj - 1
-                           zw         = zAV(ji,jj) / zBV(ji,jj-1)
+                           zfac             = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) )
+                           zw               = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 )
                            zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1)
                            zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 
@@ -821,9 +1014,14 @@
             
                         ! Backward sweep
-                        v_ice(ji,jpj-1)  = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1)  ! last row
+                        zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) )
+                        v_ice(ji,jpj-1)  = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) & 
+                               &         * vmask(ji,jpj-1,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,1) / zBV_prime(ji,jj)
+                           zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) )
+                           v_ice(ji,jj)   = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) ) &
+                               &                             / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )
                         END DO            
-            
+
                         !---------------       
                         ! -- Relaxation
@@ -835,5 +1033,5 @@
                            v_ice(ji,jj) =   zmsk00y(ji,jj)                                        & 
                               &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)             &
-                                          + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    )
+                              &           + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) * vmask(ji,jj,1)
 
                         END DO
@@ -844,6 +1042,6 @@
                                     
                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. )
+
+               CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
                               
                !--------------------------------------------------------------------------------------
@@ -856,4 +1054,5 @@
                ! 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                            
+               zuerr_max = 0._wp
                                
                IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
@@ -886,4 +1085,5 @@
                ! on V
                !------
+               zverr_max = 0._wp
                
                IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
@@ -919,19 +1119,64 @@
                !
                !---------------------------------------------------------------------------------------
-               IF( ln_rhg_chkcvg ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, &
+               IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) 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 )
-               
-            
+
+               IF ( lwp ) WRITE(numout,*) ' Done convergence tests '
+               
             ENDIF ! ---    end ll_u_iterate or ll_v_iterate
-            
+
          END DO ! i_inn, end of inner loop
 
+      END DO ! End of outer loop (i_out) =============================================================================================
+
+      IF ( lwp ) WRITE(numout,*) ' We are out of outer loop '
+
+      CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG
+      CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG
+      CALL iom_put( 'zFU_prime'     , zFU_prime      ) ! MV DEBUG
+
+      CALL iom_put( 'zFV'           , zFV            ) ! MV DEBUG
+      CALL iom_put( 'zBV_prime'     , zBV_prime      ) ! MV DEBUG
+      CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG
+
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
+
+      IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg '
+      IF( iom_use('uice_dbg' ) )   CALL iom_put( 'uice_dbg'   , u_ice    )                              ! ice velocity u after solver
+      IF( iom_use('vice_dbg' ) )   CALL iom_put( 'vice_dbg'   , v_ice    )                              ! ice velocity v after solver
+            
       !------------------------------------------------------------------------------!
       !
-      ! 6) Mask final velocities
+      ! --- Convergence diagnostics 
       !
       !------------------------------------------------------------------------------!
 
-      END DO ! End of outer loop (i_out) =============================================================================================
+      IF( ln_rhg_chkcvg ) THEN
+          
+         IF( iom_use('uice_cvg')  ) THEN
+            CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , &                ! ice velocity difference at last iteration
+                  &                        ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
+         ENDIF   
+        
+      ENDIF ! ln_rhg_chkcvg
+
+      ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead
+      DO jj = 2, jpj - 1
+         DO ji = 2, jpi - 1   
+
+             u_ice(ji,jj) =   zmsk00x(ji,jj)                               & 
+      &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    &
+                  + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    )
+
+             v_ice(ji,jj) =   zmsk00y(ji,jj)                               & 
+      &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    &
+                  + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    )
+
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
+
+      ! END DEBUG
 
       !------------------------------------------------------------------------------!
@@ -942,5 +1187,5 @@
       !
       ! MV OPT: subroutinize ?
-            
+
       DO jj = 1, jpj - 1
          DO ji = 1, jpi - 1
@@ -1087,4 +1332,5 @@
       IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
 
@@ -1123,7 +1369,8 @@
       ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765)
       ! R2 means we need to recompute stresses
+
       IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
          !
-         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )         
+         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !         
          DO jj = 2, jpj - 1
@@ -1142,7 +1389,7 @@
 
                ! Normalized  principal stresses (used to display the ellipse)
-               z1_strength      =   1._wp / strength(ji,jj)
-               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  + zsig_II(ji,jj) ) * z1_strength
-               zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) - zsig_II(ji,jj) ) * z1_strength
+               z1_strength      =   1._wp / MAX ( 1._wp , strength(ji,jj) )
+               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
             END DO
          END DO
@@ -1154,5 +1401,5 @@
          CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
 
-         DEALLOCATE( zsig1_p , zsig2_p )
+         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
          
       ENDIF
@@ -1253,15 +1500,4 @@
       ENDIF
 
-      ! --- Convergence diagnostics --- !
-      IF( ln_rhg_chkcvg ) THEN
-          
-         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
-
-      !
       DEALLOCATE( zmsk00, zmsk15 )
 
@@ -1285,5 +1521,4 @@
       !!                   - 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)   
@@ -1296,7 +1531,8 @@
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV
       !!
-      INTEGER           ::   it, idtime, istatus
+      INTEGER           ::   it, idtime, istatus, ix_dim, iy_dim
       INTEGER           ::   ji, jj          ! dummy loop indices
       REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars
+      REAL(wp)          ::   z1_pglob_area
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays
                                                                              
@@ -1304,19 +1540,27 @@
       !!----------------------------------------------------------------------
 
+      IF( lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) ' kiter    =  : ', kiter
+         WRITE(numout,*) ' kitermax =  : ', kitermax
+      ENDIF
+
       ! 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_DIM( ncvgid, 'x'     , jpi, ix_dim )
+            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim )
+
+            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid_ucvg ) ! the name uice_cvg sucks, no ?
+            ! i suggest vel_dif instead
             istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE   , (/ idtime /), nvarid_ures )
             istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE   , (/ idtime /), nvarid_vres )
@@ -1325,30 +1569,32 @@
             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_DEF_VAR( ncvgid, 'usub'    , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_usub) ! --> u-velocity in sub-iterations
+            istatus = NF90_DEF_VAR( ncvgid, 'vsub'    , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_vsub) ! --> v-velocity in sub-iterations
+
             istatus = NF90_ENDDEF(ncvgid)
+
          ENDIF
          !
       ENDIF
 
-      ! time
-      it = ( kt - 1 ) * kitermax + kiter
+      IF ( lwp ) WRITE(numout,*) ' File created '
 
       ! --- 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)
-      
+
+      IF ( lwp ) WRITE(numout,*) ' zfeldif : ', zveldif
+      
+      ! ---  Mean residual and kinetic energy
+      IF ( kiter == 1 ) THEN
+
+              zu_res_mean = 0._wp
+              zv_res_mean = 0._wp
+              zvelres     = 0._wp
+              zmke        = 0._wp
+
+      ELSE
+
       ! -- Mean residual (N/m^2), zu_res_mean
       ! Here we take the residual of the linear system (N/m^2), 
@@ -1356,18 +1602,30 @@
       ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 
       ! i.e., how close we are to a solution
+
+      IF ( lwp ) WRITE(numout,*) ' TEST 1 ' 
+
+      z1_pglob_area = 1._wp / pglob_area
+
       DO jj = 2, jpj - 1
          DO ji = 2, jpi - 1                                      
-            zu_res(ji,jj)  = prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) &
-               &           - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj)
+            zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               &
+               &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
                            
-            zv_res(ji,jj)  = prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) &
-               &           - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1)                           
+            zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               &
+               &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
+
+            zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area
+            zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area
+
           END DO
        END DO                  
-       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:,1) )
-       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:,1) )
-       zu_res_mean = SQRT( zu_res_mean / pglob_area )
-       zv_res_mean = SQRT( zv_res_mean / pglob_area )
+
+       IF ( lwp ) WRITE(numout,*) ' TEST 2 ' 
+       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) )
+       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) )
+       IF ( lwp ) WRITE(numout,*) ' TEST 3 ' 
        zvelres     = 0.5_wp * ( zu_res_mean + zv_res_mean )
+
+       IF ( lwp ) WRITE(numout,*) ' TEST 4 ' 
                                           
        ! -- Global mean kinetic energy per unit area (J/m2)  
@@ -1380,17 +1638,29 @@
        END DO
        
-       zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) )
+         IF ( lwp ) WRITE(numout,*) ' TEST 5 ' 
+       zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area
+         IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 
+
+      ENDIF
                   
          !                                                ! ==================== !
 
+      ! time
+      it = ( kt - 1 ) * kitermax + kiter
+
+
       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/) )
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) )     ! max U or V velocity diff between subiterations
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) )   ! average of u- and v- residuals
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) )   ! max U velocity difference, inner iterations
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) )   ! max V velocity difference, inner iterations
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) )         ! mean kinetic energy
+
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_usub, (/pu/), (/it/) )          ! u-velocity
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_vsub, (/pv/), (/it/) )          ! v-velocity
+
          ! close file
          IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid )
