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 15526)
+++ /NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 15527)
@@ -152,5 +152,5 @@
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass and volume
       REAL(wp) ::   zds2, zdt, zdt2, zdiv, zdiv2                        ! temporary scalars
-      REAL(wp) ::   zp_deltastar_f                                      ! 
+      REAL(wp) ::   zp_delstar_f                                      ! 
       REAL(wp) ::   zu_cV, zv_cU                                        ! 
       REAL(wp) ::   zfac, zfac1, zfac2, zfac3
@@ -164,7 +164,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat, zdeltastar_t           ! Delta & Delta* at T-points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat, zdelstar_t           ! Delta & Delta* at T-points
       REAL(wp), DIMENSION(jpi,jpj) ::   ztens, zshear                   ! Tension, shear
-      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta* at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delstar_t                  ! P/delta* at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zef                             ! Viscosity pre-factor at F point
@@ -194,5 +194,4 @@
       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)
@@ -302,6 +301,4 @@
          END DO
       END DO
-
-      !CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp ) ! OK to remove
 
       ! Lateral boundary conditions on velocity (modify zfmask)
@@ -434,5 +431,5 @@
          END DO
 
-         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! necessary 
+         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! necessary, zds2 uses jpi/jpj values for zds 
 
          DO jj = 2, jpj        ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
@@ -460,11 +457,11 @@
                
                ! delta* at T points (following Lemieux and Dupont, GMD 2020)
-               zdeltastar_t(ji,jj)   = zdeltat(ji,jj) + rn_creepl
+               zdelstar_t(ji,jj)     = zdeltat(ji,jj) + rn_creepl ! OPT zdelstar_t can be totally removed and put into next line directly. Could change results
               
-               ! P/delta at T-points
-               zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj)
+               ! P/delta* at T-points
+               zp_delstar_t(ji,jj)   = strength(ji,jj) / zdelstar_t(ji,jj)
                
                ! Temporary zzt and zet factors at T-points
-               zzt(ji,jj)            = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)
+               zzt(ji,jj)            = zp_delstar_t(ji,jj) * r1_e1e2t(ji,jj)
                zet(ji,jj)            = zzt(ji,jj)     * z1_ecc2 
                           
@@ -472,7 +469,5 @@
          END DO
          
-!         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. , zdeltat, 'T', 1.) ! OK to
-!         remove zet and zzt
-         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1.,  zdeltat, 'T', 1.)
+         CALL lbc_lnk( 'icedyn_rhg_vp', zp_delstar_t , 'T', 1. ) ! necessary, used for ji = 1 and jj = 1
 
          DO jj = 1, jpj - 1
@@ -480,13 +475,11 @@
                
                ! 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) )
+               zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) )
                
                ! Temporary zef factor at F-point
-               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) * 0.5_wp ! BUG missing un facteur 2 Nov 12
-
-            END DO
-         END DO
-
-!        CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) ! Ok to remove
+               zef(ji,jj)      = zp_delstar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) * 0.5_wp 
+
+            END DO
+         END DO
 
          !---------------------------------------------------
@@ -524,10 +517,4 @@
          END DO
 
-!        CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCwU ,  'U', -1., zCwV, 'V', -1. )! OK to remove
-!        CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCorU,  'U', -1., zCorV, 'V', -1. )! OK to remove
-
-         ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system, 
-         ! so there is no need for lbclnk on drag and coriolis
-
          !-------------------------------------
          ! -- Internal stress RHS contribution
@@ -540,5 +527,5 @@
                ! sig1 contribution to RHS of U-equation at T-points
                zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_c(ji,jj) - e1v(ji,jj-1)    * zv_c(ji,jj-1) )   &
-                               &                - zp_deltastar_t(ji,jj) * zdeltat(ji,jj)
+                               &                - zp_delstar_t(ji,jj) * zdeltat(ji,jj)
                                             
                ! sig2 contribution to RHS of U-equation at T-points            
@@ -547,5 +534,5 @@
                ! sig1 contribution to RHS of V-equation at T-points
                zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj)    * zu_c(ji,jj) - e2u(ji-1,jj)    * zu_c(ji-1,jj) )   & 
-                               &                - zp_deltastar_t(ji,jj) * zdeltat(ji,jj)
+                               &                - zp_delstar_t(ji,jj) * zdeltat(ji,jj)
 
                ! sig2 contribution to RHS of V-equation  at T-points
@@ -554,6 +541,4 @@
             END DO
          END DO
-
-         ! a priori, no lbclnk, because rhsu is only used in the inner domain
          
          ! --- Stress contributions at F-points         
@@ -572,9 +557,4 @@
             END DO
          END DO
-
-         ! CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) ! Ok to remove
-         ! CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) ! Ok to remove
-
-         ! a priori, no lbclnk, because rhsu are only used in the inner domain
 
          ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002)
@@ -613,8 +593,4 @@
          END DO
          
-!        CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V',  -1.)! OK to remove
-!        CALL lbc_lnk_multi( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V',  -1.)! OK to remove
-!        CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V',  -1.)! OK to remove
-
          !---------------------------------------------------------------------------------------!
          !
@@ -707,11 +683,4 @@
          END DO
 
-         ! OK to remove
-         !CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU  , 'U', 1., zAV  , 'V',  1. ) ! --> here normal that we use '1' and not '-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. )
-
       !------------------------------------------------------------------------------!
       !
@@ -745,5 +714,5 @@
                   ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F
 
-                  zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; zCU_prime(:,:) = 0._wp
+                  zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; 
                   
                   DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! )
@@ -761,6 +730,8 @@
                         !------------------------
                         DO ji = 2, jpi - 1    
-
-                           ! boundary condition substitution
+                           ! note: these are key lines linking information between processors
+                           ! u_ice/v_ice need to be lbc_linked
+
+                           ! sub-domain boundary condition substitution
                            ! see Zhang and Hibler, 1997, Appendix B
                            zAA3 = 0._wp
@@ -778,6 +749,4 @@
                      END DO
                      
-                     CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U',  -1. ) ! -1, zFU is a vector
-                     
                      !---------------
                      ! Forward sweep
@@ -785,5 +754,4 @@
                      DO jj = jj_min, jpj - 1, nn_zebra_vp
       
-                        ! MV BUG OCT 6. Both zBU_prime and zFU_prime were undefined for ji = 2
                         zBU_prime(2,jj)     = zBU(2,jj)
                         zFU_prime(2,jj)     = zFU(2,jj)
@@ -799,5 +767,4 @@
 
                      END DO
-                     ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime, 'U',  -1., zBU_prime, 'U',  1. ) --> can be removed
                                                                                                      
                      !-----------------------------
@@ -820,6 +787,5 @@
                         END DO
 
-                        !--- Relaxation
-                        ! and velocity masking for little-ice and no-ice cases
+                        !--- Relaxation and masking (for low-ice/no-ice cases)
                         DO ji = 2, jpi - 1    
                         
@@ -848,5 +814,5 @@
                   !!! ZH97 explain it is critical for convergence speed
 
-                  zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp
+                  zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; 
 
                   DO jn = 1, nn_zebra_vp ! "zebra" loop
@@ -863,5 +829,5 @@
                         DO jj = 2, jpj - 1
 
-                           ! boundary condition substitution (check it is correctly applied !!!)
+                           ! subdomain boundary condition substitution (check it is correctly applied !!!)
                            ! see Zhang and Hibler, 1997, Appendix B
                            zAA3 = 0._wp
@@ -879,7 +845,4 @@
                      END DO
 
-                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV, 'V',  -1.) ! MV Nov 5 bug test to change boundary condition to -1
-                                                                          ! since zFV is a vector!!!
-                     
                      !---------------
                      ! Forward sweep
@@ -887,5 +850,4 @@
                      DO ji = ji_min, jpi - 1, nn_zebra_vp 
                      
-                        ! MV BUG OCT 6.. Both zBV_prime and zFV_prime were undefined for jj = 2
                         zBV_prime(ji,2)     = zBV(ji,2)
                         zFV_prime(ji,2)     = zFV(ji,2)
@@ -902,7 +864,4 @@
                      END DO
 
-                     ! new test
-                     ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV_prime, 'V',  -1., zBV_prime, 'V',  1. ) ! MV Nov 5 bug test
-                     
                      !-----------------------------
                      ! Backward sweep & relaxation
@@ -923,5 +882,5 @@
                         END DO            
                                                    
-                        ! --- Relaxation  & masking (should it be now or later)
+                        ! --- Relaxation & masking 
                         DO jj = 2, jpj - 1
                         
@@ -940,4 +899,5 @@
                ENDIF !   ll_v_iterate
 
+               ! I suspect the communication should go into the zebra loop if we want reproducibility
                CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
                               
@@ -1031,6 +991,4 @@
 
       END DO ! End of outer loop (i_out) =============================================================================================
-
-      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
 
       IF( ln_rhg_chkcvg ) THEN
@@ -1123,6 +1081,6 @@
          DO jj = 2, jpj - 1
             DO ji = 2, jpi - 1
-               zp_deltastar_t(ji,jj)   =   strength(ji,jj) / pdelta_i(ji,jj) 
-               zfac                    =   zp_deltastar_t(ji,jj) 
+               zp_delstar_t(ji,jj)     =   strength(ji,jj) / pdelta_i(ji,jj) 
+               zfac                    =   zp_delstar_t(ji,jj) 
                zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
                zs2(ji,jj)              =   zfac * z1_ecc2 * ztens(ji,jj)
@@ -1142,8 +1100,8 @@
             
                ! 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) )
+               zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) )
                
                ! s12 at F-points 
-               zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj)
+               zs12f(ji,jj) = zp_delstar_f * z1_ecc2 * zds(ji,jj)
                
             END DO
@@ -1251,5 +1209,5 @@
                !                        and **deformations** at current iterates
                !                        following Lemieux & Dupont (2020)
-               zfac             =   zp_deltastar_t(ji,jj)
+               zfac             =   zp_delstar_t(ji,jj)
                zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) )
                zsig2            =   zfac * z1_ecc2 * ztens(ji,jj)
@@ -1618,5 +1576,4 @@
          ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 
          ! i.e., how close we are to a solution
-         !!!! UNITS OF RESIDUAL ARE NOT m/s BUT N/m2 (CHECK)
          zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp
 
