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 13591)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13593)
@@ -131,5 +131,5 @@
       LOGICAL ::   ll_u_iterate, ll_v_iterate   ! continue iteration or not
       !
-      INTEGER ::   ji, ji2, jj, jn          ! dummy loop indices
+      INTEGER ::   ji, ji2, jj, jj2, jn          ! dummy loop indices
       INTEGER ::   jter, i_out, i_inn  ! 
       INTEGER ::   ji_min, jj_min      !
@@ -777,11 +777,6 @@
          ll_u_iterate = .TRUE.
          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
+         
+         nn_thomas = 2 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm
 
          DO i_inn = 1, nn_ninn_vp ! inner loop iterations
@@ -806,6 +801,8 @@
                                            ! ---------------------------- !
       
+                  ! What follows could be subroutinized...
+      
                   ! Thomas Algorithm  for tridiagonal solver
-                  ! what follows could be subroutinized...
+                  ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F
                   
                   DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! )
@@ -819,21 +816,13 @@
                      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
-
-                        !
-                        ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan
-                        !
-                        ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F
-
-                        ! - Right-hand side of tridiagonal system (zFU)
+                     
+                        !------------------------
+                        ! Independent term (zFU)
+                        !------------------------
                         DO ji = 2, jpi - 1    
 
                            ! boundary condition substitution
                            ! see Zhang and Hibler, 1997, Appendix B
-                           ! MV NOTE possibly not fully appropriate
                            zAA3 = 0._wp
                            IF ( ji == 2 )          zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj)
@@ -848,16 +837,9 @@
                         END DO
 
-                     END DO
-
-                     ! CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. )
-                        
-                     !---------------
-                     ! 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
+                        !---------------
+                        ! Forward sweep
+                        !---------------
+         
+                        IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm
 
                            DO ji = 3, jpi - 1
@@ -870,9 +852,5 @@
                            END DO
                         
-                        END DO
-
-                     ELSE ! nn_thomas == 2      ! mitGCM algorithm
-
-                        DO jj = jj_min, jpj - 1, nn_zebra_vp
+                        ELSE ! nn_thomas == 2      ! mitGCM algorithm
 
                            ! ji = 2
@@ -893,16 +871,12 @@
                            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
+                        ENDIF
+ 
+                        !-----------------------------
+                        ! Backward sweep & relaxation
+                        !-----------------------------
+                    
+                        ! --- Backward sweep 
+                        IF ( nn_thomas == 1 ) THEN ! MV version
 
                            ! last row 
@@ -912,5 +886,5 @@
 
                            DO ji2 = 2 , jpi-2 ! all other rows       
-                          !DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop
+                           !DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop
                               ji =  jpi - ji2
                               ! ji2 = 2       -> ji = jpi - 2
@@ -922,11 +896,7 @@
                            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)
+                        ELSE ! nn_thomas == 2 ! mitGCM version
+
+                           u_ice(jpi-1,jj)   =   zFU_prime(jpi-1,jj)
 
                            DO ji2 = 2 , jpi-2 ! all other rows       !   MV OPT: could be turned into forward loop (by substituting ji)
@@ -935,20 +905,13 @@
                            END DO
                      
-                        END DO
-
-                     ENDIF
-
-                     !---------------       
-                     ! -- Relaxation
-                     !---------------
-
-                     DO jj = jj_min, jpj - 1, nn_zebra_vp
-            
+                        ENDIF
+
+                        !--- Relaxation
+                        ! and velocity masking for little-ice and no-ice cases
                         DO ji = 2, jpi - 1    
                         
-                           u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) )
+                           u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) ! relaxation
                            
-                           ! velocity masking for little-ice and no-ice cases
-                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        & 
+                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        &   ! masking
                               &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             &
                               &           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) * umask(ji,jj,1)
@@ -966,6 +929,10 @@
                !                           ! ---------------------------- !
                                           
-                  ! MV OPT: what follows could be subroutinized...
-                  
+                  ! MV OPT: what follows could be subroutinized...                  
+                  ! Thomas Algorithm  for tridiagonal solver
+                  ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F
+                  ! It is intentional to have a ji then jj loop for V-velocity
+                  !!! ZH97 explain it is critical for convergence speed
+
                   DO jn = 1, nn_zebra_vp ! "zebra" loop
       
@@ -977,14 +944,8 @@
          
                      DO ji = ji_min, jpi - 1, nn_zebra_vp 
-                        
-                        !!! It is intentional to have a ji then jj loop for V-velocity
-                        !!! ZH97 explain it is critical for convergence speed
-
-                        !-------------------------------------------
-                        ! -- Tridiagonal system solver for each row
-                        !-------------------------------------------
-                        ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F
-
-                        ! --- Right-hand side of tridiagonal system (zFU)
+                     
+                        !------------------------
+                        ! Independent term (zFU)
+                        !------------------------
                         DO jj = 2, jpj - 1
 
@@ -1002,41 +963,86 @@
 
                         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
-                           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
-            
-                        ! Backward sweep
-                        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
-                           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
                         !---------------
+                     
+                        IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than 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
+																											      
+                        ELSE ! nn_thomas == 2      ! mitGCM algorithm
+                     
+                           ! jj = 2
+                           zw                   =   zBV(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 )
+                           zCV_prime(2,jj)      =   zw * zCV(2,jj)
+                           zFV_prime(2,jj)      =   zw * zFV(2,jj)
+
+                           DO jj = 3, jpj - 1
+
+                              zw                =   zBV(ji,jj) - zAV(ji,jj) * zCV_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 )
+                              zCV_prime(ji,jj)  =   zw * zCV(ji,jj) 
+                              zFV_prime(ji,jj)  =   zw * ( zFV(ji,jj) - zAV(ji,jj) * zFV_prime(ji-1,jj) )
+
+                           END DO
+                       
+                        ENDIF
+                     
+                        !-----------------------------
+                        ! Backward sweep & relaxation
+                        !-----------------------------
+                    
+                        ! --- 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)
+                           
+                           ! other rows
+                           DO jj2 = 2 , jpi-2 ! all other rows
+                              jj =  jpj - jj2
+                              v_ice(ji,jj)   =   zFV_prime(ji,jj) - zCV_prime(ji,jj) * v_ice(ji,jj+1)
+                           END DO
+                        
+                        ENDIF ! nn_thomas
+                        
+                        ! --- Relaxation  & masking (should it be now or later)
                         DO jj = 2, jpj - 1
-                           v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) )
-
-                           ! mask velocity for no-ice and little-ice cases                           
-                           v_ice(ji,jj) =   zmsk00y(ji,jj)                                        & 
-                              &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)             &
+                        
+                            v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) )    ! relaxation
+                            
+                            v_ice(ji,jj) =   zmsk00y(ji,jj)                                        &      ! masking
+                              &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)              &
                               &           + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) * vmask(ji,jj,1)
 
-                        END DO
-
+                        END DO ! jj
+                        
                      END DO ! ji
-
+                     
                   END DO ! zebra loop
                                     
