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 13627)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13629)
@@ -777,6 +777,7 @@
          ll_u_iterate = .TRUE.
          ll_v_iterate = .TRUE.
-         
-         nn_thomas = 2 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm
+
+!        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
@@ -789,9 +790,9 @@
             ! 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
+            zu_b(:,:)      = u_ice(:,:) ! velocity at previous sub-iterate
+            zv_b(:,:)      = v_ice(:,:)
+
+!           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
 
             IF ( ll_u_iterate .OR. ll_v_iterate )   THEN
@@ -800,9 +801,11 @@
                IF ( ll_u_iterate ) THEN    ! --- Solve for u-velocity --- !
                                            ! ---------------------------- !
-      
+
                   ! What follows could be subroutinized...
       
                   ! Thomas Algorithm  for tridiagonal solver
                   ! 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
                   
                   DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! )
@@ -826,6 +829,6 @@
                            ! see Zhang and Hibler, 1997, Appendix B
                            zAA3 = 0._wp
-                           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)
+                           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
@@ -854,16 +857,17 @@
                         ELSE ! nn_thomas == 2      ! mitGCM algorithm
 
-                           ! ji = 2
+                           ! --- 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 )
+                           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 * 1.0_wp /  MAX ( 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) )
@@ -882,5 +886,5 @@
                            ! 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 ) & 
+                           u_ice(jpi-1,jj)    = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 
                                   &            * umask(jpi-1,jj,1)
 
@@ -891,5 +895,5 @@
                               ! 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
+                              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 ) 
@@ -898,9 +902,11 @@
                         ELSE ! nn_thomas == 2 ! mitGCM version
 
-                           u_ice(jpi-1,jj)   =   zFU_prime(jpi-1,jj)
+                           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
-                              u_ice(ji,jj)    =   zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj)
+                              ! 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
                      
@@ -915,5 +921,5 @@
                            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)
+                              &           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp     ) * umask(ji,jj,1)
                            
                         END DO
@@ -935,4 +941,6 @@
                   !!! ZH97 explain it is critical for convergence speed
 
+                  zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp
+
                   DO jn = 1, nn_zebra_vp ! "zebra" loop
       
@@ -946,5 +954,5 @@
                      
                         !------------------------
-                        ! Independent term (zFU)
+                        ! Independent term (zFV)
                         !------------------------
                         DO jj = 2, jpj - 1
@@ -970,30 +978,28 @@
                         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
-																											      
+                           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)
+                           zw                   =   zBV(ji,2)
                            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)
+                           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-1,jj)
+                              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 * 1.0_wp /  MAX ( 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-1,jj) )
-
+                              zFV_prime(ji,jj)  =   zw * ( zFV(ji,jj) - zAV(ji,jj) * zFV_prime(ji,jj-1) )
                            END DO
-                       
+
                         ENDIF
                      
@@ -1015,5 +1021,5 @@
                               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) ) &
+                              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            
@@ -1022,12 +1028,12 @@
 
                            ! last row                        
-                           v_ice(ji,jpj-1)   =   zFV_prime(ji,jpj)
+                           v_ice(ji,jpj-1)   =   zFV_prime(ji,jpj-1) * vmask(ji,jpj-1,1)
                            
                            ! other rows
-                           DO jj2 = 2 , jpi-2 ! all 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)
+                              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
                         
@@ -1138,10 +1144,17 @@
       IF ( lwp ) WRITE(numout,*) ' We are out of outer loop '
 
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU  , 'U',  1., zFV  , 'V',  1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU_prime  , 'U',  1., zBV_prime  , 'V',  1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime  , 'U',  1., zFV_prime  , 'V',  1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU_prime  , 'U',  1., zCV_prime  , 'V',  1. )
+
       CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG
       CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG
+      CALL iom_put( 'zCU_prime'     , zCU_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( 'zCV_prime'     , zCV_prime      ) ! MV DEBUG
       CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG
 
@@ -1183,4 +1196,6 @@
 
       CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
+
+      IF ( lwp ) WRITE(numout,*) ' Velocity replaced '
 
       ! END DEBUG
@@ -1238,4 +1253,6 @@
          END DO
       END DO
+
+      IF ( lwp ) WRITE(numout,*) ' Deformation recalculated '
       
       CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
@@ -1288,4 +1305,7 @@
          
       ENDIF
+
+      IF ( lwp ) WRITE(numout,*) ' zs12f recalculated '
+
       !
       !-----------------------
@@ -1341,4 +1361,6 @@
       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
 
+      IF ( lwp ) WRITE(numout,*) 'Some terms recalculated '
+
       ! --- Stress tensor invariants (SIMIP diags) --- !
       IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
@@ -1368,4 +1390,6 @@
          
       ENDIF
+
+      IF ( lwp ) WRITE(numout,*) 'SIMIP work done'
 
       ! --- Normalized stress tensor principal components --- !
@@ -1387,4 +1411,5 @@
                zfac             =   zp_deltastar_t(ji,jj)
                zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) )
+               zsig1            = 0._wp !!! FUCKING DEBUG TEST !!!
                zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
                zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
@@ -1400,7 +1425,9 @@
             END DO
          END DO
+         IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done'
          !
          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
          !      
+         IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll '
          !
          CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
@@ -1408,4 +1435,6 @@
 
          DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
+
+         IF ( lwp ) WRITE(numout,*) ' So what ??? '
          
       ENDIF
