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 13680)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13681)
@@ -43,5 +43,4 @@
    !! for convergence tests
    INTEGER ::   ncvgid        ! netcdf file id
-   INTEGER ::   nvarid_ucvg   ! netcdf variable id
    INTEGER ::   nvarid_ures
    INTEGER ::   nvarid_vres
@@ -49,6 +48,7 @@
    INTEGER ::   nvarid_udif
    INTEGER ::   nvarid_vdif
+   INTEGER ::   nvarid_veldif
    INTEGER ::   nvarid_mke
-   INTEGER ::   nvarid_usub, nvarid_vsub
+   INTEGER ::   nvarid_ures_xy, nvarid_vres_xy
 
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
@@ -136,6 +136,4 @@
       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)      
       !
@@ -144,5 +142,5 @@
       REAL(wp) ::   zglob_area                                          ! global ice area for diagnostics
       REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
-      REAL(wp) ::   zm2, zm3, zmassU, zmassV                            ! ice/snow mass and volume
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass and volume
       REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars
       REAL(wp) ::   zp_deltastar_f                                      ! 
@@ -294,5 +292,4 @@
          END DO
       END DO
-      CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp )
 
       ! Lateral boundary conditions on velocity (modify zfmask)
@@ -335,4 +332,7 @@
       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
       zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!!
+
+      zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:)       ! Snow and ice mass at T-point
+      zmf(:,:) = zmt(:,:) * ff_t(:,:)                      ! Coriolis factor at T points (m*f)
       
       DO jj = 2, jpj - 1
@@ -344,10 +344,10 @@
 
             ! Snow and ice mass at U-V points
-            zmt(ji,jj)      = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
-            zm2             = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
-            zm3             = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
+            zm1             = zmt(ji,jj)
+            zm2             = zmt(ji+1,jj)
+            zm3             = zmt(ji,jj+1)
             
-            zmassU          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
-            zmassV          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+            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 area divided by time step
@@ -362,7 +362,4 @@
             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)      = zmt(ji,jj) * ff_t(ji,jj)
             
             ! Wind stress
@@ -379,8 +376,18 @@
 
             ! 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
-            IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
-            ELSE                                                     ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
+            IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+            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              
+
+! MV TEST DEBUG
+            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji+1,jj)  <= zmmin )     .AND.  &
+               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji+1,jj) <= zamin ) )    THEN   ;   zmsk01x(ji,jj) = 0._wp
+            ELSE                                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+
+            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji,jj+1)  <= zmmin )     .AND.  &
+               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji,jj+1) <= zamin ) )    THEN   ;   zmsk01y(ji,jj) = 0._wp
+            ELSE                                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
+! END MV TEST DEBUG
 
          END DO
@@ -391,5 +398,4 @@
       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
@@ -449,6 +455,11 @@
          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
-            DO ji = 2, jpi - 1 ! 
+         !DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
+         !   DO ji = 2, jpi - 1 ! 
+
+! MV DEBUG
+         DO jj = 2, jpj        ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
+            DO ji = 2, jpi     ! 
+! END MV DEBUG
 
                ! shear**2 at T points (doc eq. A16)
@@ -500,5 +511,5 @@
                
                ! Temporary zef factor at F-point
-               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2
+               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj)
 
             END DO
@@ -545,12 +556,13 @@
 
          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 lbc_lnk_multi( 'icedyn_rhg_vp', zCwU ,  'U', -1., zCwV, 'V', -1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCorU,  'U', -1., zCorV, '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
          
@@ -652,6 +664,12 @@
          END DO
          
-         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', 1., zrhsv, 'V',  1. )
-
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V',  -1.)
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V',  -1.)
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V',  -1.)
+
+         CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG
+         CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG
+         CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG
+         CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG
          CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG
          CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG
@@ -778,7 +796,4 @@
          ll_v_iterate = .TRUE.
 
-!        nn_thomas = 2 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm ! --->  not working for now
-         nn_thomas = 1 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm ! ---> BUG MPP!!!
-
          DO i_inn = 1, nn_ninn_vp ! inner loop iterations
 
@@ -840,76 +855,42 @@
                         END DO
 
-                        !---------------
-                        ! Forward sweep
-                        !---------------
-         
-                        IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm
-
-                           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
-                        
-                        ELSE ! nn_thomas == 2      ! mitGCM algorithm
-
-                           ! --- ji = 2
-                           zw                   =   zBU(2,jj)
-                           zfac                 =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) )
-                           zw                   =   zfac /  MAX ( ABS( zw ) , epsi20 )
-                           zCU_prime(2,jj)      =   zw * zCU(2,jj)
-                           zFU_prime(2,jj)      =   zw * zFU(2,jj)
-
-                           ! --- ji = 3 --> jpi-1
-                           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 /  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
-
-                        ENDIF
+                     END DO
+                     
+                     CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U',  1. )
+                     
+                     !---------------
+                     ! Forward sweep
+                     !---------------   
+                     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
+
+                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime, 'U',  1., zBU_prime, 'U', 1. )
  
-                        !-----------------------------
-                        ! Backward sweep & relaxation
-                        !-----------------------------
+                     !-----------------------------
+                     ! Backward sweep & relaxation
+                     !-----------------------------
+
+                     DO jj = jj_min, jpj - 1, nn_zebra_vp
                     
                         ! --- Backward sweep 
-                        IF ( nn_thomas == 1 ) THEN ! MV version
-
-                           ! 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            
-
-                        ELSE ! nn_thomas == 2 ! mitGCM version
-
-                           u_ice(jpi-1,jj)    =   zFU_prime(jpi-1,jj) * umask(jpi-1,jj,1)
-
-                           DO ji2 = 2 , jpi-2 ! all other rows       !   MV OPT: could be turned into forward loop (by substituting ji)
-                              ji =  jpi - ji2
-                              ! ji2 = 2, ji = jpi - 2
-                              ! ji2 = jpi- 2, ji = 2
-                              u_ice(ji,jj)    = ( zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1)
-                           END DO
-                     
-                        ENDIF
+                        ! 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 ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop
+                           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)   & 
+                                           &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) 
+                        END DO
 
                         !--- Relaxation
@@ -972,70 +953,44 @@
                         END DO
                         
-                        !---------------
-                        ! Forward sweep
-                        !---------------
+                     END DO
+
+                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV, 'V',  1.)
                      
-                        IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm
+                     !---------------
+                     ! Forward sweep
+                     !---------------
+                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
                      
-                           DO jj = 3, jpj - 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) 
-                           END DO
-
-                        ELSE ! nn_thomas == 2      ! mitGCM algorithm
+                        DO jj = 3, jpj - 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) 
+
+                        END DO
+
+                     END DO
+
+                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV_prime, 'V',  1., zBV_prime, 'V', 1. )
                      
-                           ! jj = 2
-                           zw                   =   zBV(ji,2)
-                           zfac                 =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) )
-                           zw                   =   zfac /  MAX ( ABS( zw ) , epsi20 )
-                           zCV_prime(ji,2)      =   zw * zCV(ji,2)
-                           zFV_prime(ji,2)      =   zw * zFV(ji,2)
-
-                           DO jj = 3, jpj - 1
-                              zw                =   zBV(ji,jj) - zAV(ji,jj) * zCV_prime(ji,jj-1)
-                              zfac              =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) )
-                              zw                =   zfac /  MAX ( ABS( zw ) , epsi20 )
-                              zCV_prime(ji,jj)  =   zw * zCV(ji,jj) 
-                              zFV_prime(ji,jj)  =   zw * ( zFV(ji,jj) - zAV(ji,jj) * zFV_prime(ji,jj-1) )
-                           END DO
-
-                        ENDIF
-                     
-                        !-----------------------------
-                        ! Backward sweep & relaxation
-                        !-----------------------------
+                     !-----------------------------
+                     ! Backward sweep & relaxation
+                     !-----------------------------
+                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
                     
                         ! --- Backward sweep 
-                        IF ( nn_thomas == 1 ) THEN ! nn_thomas = 1    ! --- MV version seemingly more than mitGCM algorithm
-             
-                           ! 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
-
-                           ! other rows
-                           !DO jj = jpj-2, 2, -1 ! original back loop
-                           DO jj2 = 2 , jpj-2    ! modified forward loop
-                              jj =  jpj - jj2    ! index swap
-                              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            
-                           
-                        ELSE                       ! nn_thomas == 2   ! --- mitGCM algorithm
-
-                           ! last row                        
-                           v_ice(ji,jpj-1)   =   zFV_prime(ji,jpj-1) * vmask(ji,jpj-1,1)
-                           
-                           ! other rows
-                           DO jj2 = 2 , jpj-2 ! all other rows
-                              jj =  jpj - jj2
-                              v_ice(ji,jj)   = ( zFV_prime(ji,jj) - zCV_prime(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1)
-                           END DO
-
-                        ENDIF ! nn_thomas
-                        
+                        ! 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
+
+                        ! other rows
+                        DO jj = jpj-2, 2, -1 ! original back loop
+                           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  & masking (should it be now or later)
                         DO jj = 2, jpj - 1
@@ -1125,16 +1080,17 @@
                   
                ENDIF ! ll_v_iterate
-               
-               !---------------------------------------------------------------------------------------
-               !
-               ! --- Calculate extra convergence diagnostics and save them
-               !
-               !---------------------------------------------------------------------------------------
-               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
+               
+            !---------------------------------------------------------------------------------------
+            !
+            ! --- Calculate extra convergence diagnostics and save them
+            !
+            !---------------------------------------------------------------------------------------
+
+            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 '
 
          END DO ! i_inn, end of inner loop
@@ -1596,15 +1552,15 @@
             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 )
-            istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE   , (/ idtime /), nvarid_velres )
-            istatus = NF90_DEF_VAR( ncvgid, 'u_dif', NF90_DOUBLE   , (/ idtime /), nvarid_udif )
-            istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE   , (/ idtime /), nvarid_vdif )
-            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE   , (/ idtime /), nvarid_mke )
-
-            istatus = NF90_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_DEF_VAR( ncvgid, 'u_res'  , NF90_DOUBLE  , (/ idtime /), nvarid_ures )
+            istatus = NF90_DEF_VAR( ncvgid, 'v_res'  , NF90_DOUBLE  , (/ idtime /), nvarid_vres )
+            istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE  , (/ idtime /), nvarid_velres )
+            istatus = NF90_DEF_VAR( ncvgid, 'u_dif'  , NF90_DOUBLE  , (/ idtime /), nvarid_udif )
+            istatus = NF90_DEF_VAR( ncvgid, 'v_dif'  , NF90_DOUBLE  , (/ idtime /), nvarid_vdif )
+            istatus = NF90_DEF_VAR( ncvgid, 'vel_dif', NF90_DOUBLE  , (/ idtime /), nvarid_veldif )
+            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE  , (/ idtime /), nvarid_mke )
+
+            istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy)
+            istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy)
 
             istatus = NF90_ENDDEF(ncvgid)
@@ -1620,13 +1576,11 @@
                                             ! 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
+         zu_res_mean = 0._wp
+         zv_res_mean = 0._wp
+         zvelres     = 0._wp
+         zmke        = 0._wp
 
       ELSE
@@ -1641,4 +1595,6 @@
 
       z1_pglob_area = 1._wp / pglob_area
+
+      zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp
 
       DO jj = 2, jpj - 1
@@ -1653,29 +1609,32 @@
             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                  
-
-       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 ' 
+         END DO
+      END DO                  
+
+      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)  
-       DO jj = 2, jpj - 1
-          DO ji = 2, jpi - 1                   
-             zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point
-             zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) )
-             zvel2(:,:)  = zu*zu + zv*zv              ! square of ice velocity at T-point  
-          END DO
-       END DO
+      ! -- Global mean kinetic energy per unit area (J/m2)
+      zvel2(:,:) = 0._wp
+      DO jj = 2, jpj - 1
+         DO ji = 2, jpi - 1                   
+            zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point
+            zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) )
+            zvel2(ji,jj)  = zu*zu + zv*zv              ! square of ice velocity at T-point  
+         END DO
+      END DO
        
-         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
+      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 ! kiter
                   
          !                                                ! ==================== !
@@ -1687,14 +1646,18 @@
       IF( lwm ) THEN
          ! write variables
-            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
+         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_veldif, (/zveldif/), (/it/), (/1/) )   ! max U or V velocity diff between subiterations
+         istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) )         ! mean kinetic energy
+
+         !
+         IF ( kiter == kitermax ) THEN
+            WRITE(numout,*) ' Should plot the spatially dependent residual '
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) )          ! U-residual, spatially dependent
+            istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) )          ! V-residual, spatially dependent
+         ENDIF
 
          ! close file
Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/iceistate.F90	(revision 13680)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/iceistate.F90	(revision 13681)
@@ -165,4 +165,6 @@
       u_ice (:,:) = 0._wp
       v_ice (:,:) = 0._wp
+
+      CALL lbc_lnk_multi( 'ice_istate' , u_ice, 'U', -1., v_ice, 'V', -1. )
       !
       !------------------------------------------------------------------------
