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 15402)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 15403)
@@ -79,5 +79,5 @@
       !!
       !!  f1(u) = g1(v)
-      !!  f2(v) = g2(v)
+      !!  f2(v) = g2(u)
       !!
       !!  The right-hand side (RHS) is explicit
@@ -143,5 +143,5 @@
       REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass and volume
-      REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars
+      REAL(wp) ::   zds2, zdt, zdt2, zdiv, zdiv2                        ! temporary scalars
       REAL(wp) ::   zp_deltastar_f                                      ! 
       REAL(wp) ::   zu_cV, zv_cU                                        ! 
@@ -156,5 +156,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat, zdeltastar_t           ! Delta & 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
@@ -207,14 +207,22 @@
       !!----------------------------------------------------------------------------------------------------------------------
       ! 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
+      ! air-ice drag
+      !utau_ice(:,:) = 0._wp
+      !vtau_ice(:,:) = 0._wp
+
+      ! coriolis
+      !ff_t(:,:) = 0._wp
+
+      ! ice-ocean drag
+      !rn_cio = 0._wp
+
+      ! sea surface height
+      !    embedded sea ice: compute representative ice top surface
+      !    non-embedded sea ice: use ocean surface for slope calculation
+      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) ! DEBUG
+      !zsshdyn(:,:) = 0._wp 
+
       ! END DEBUG
+      !!----------------------------------------------------------------------------------------------------------------------
 
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)'
@@ -265,9 +273,10 @@
       ENDIF
 
-      zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp
+      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
+      zrhsu  (:,:) = 0._wp ; zrhsv  (:,:) = 0._wp
+      zf_rhsu(:,:) = 0._wp ; zf_rhsv(:,:) = 0._wp
 
       !------------------------------------------------------------------------------!
@@ -279,7 +288,7 @@
       CALL ice_strength ! strength at T points
       
-      !------------------------------
-      ! -- F-mask       (code from EVP)
-      !------------------------------
+      !---------------------------
+      ! -- F-mask (code from EVP)
+      !---------------------------
       ! MartinV: 
       ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip)
@@ -289,5 +298,5 @@
       DO jj = 1, jpj - 1
          DO ji = 1, jpi - 1
-            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+            zfmask(ji,jj)    = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
          END DO
       END DO
@@ -326,10 +335,4 @@
       !----------------------------------------------------------------------------------------------------------
       ! Compute all terms & factors independent of velocities, or only depending on velocities at previous time step      
-      
-      ! sea surface height
-      !    embedded sea ice: compute representative ice top surface
-      !    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 !!!
 
       zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:)       ! Snow and ice mass at T-point
@@ -347,5 +350,4 @@
             zm2             = zmt(ji+1,jj)
             zm3             = zmt(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)
@@ -382,15 +384,18 @@
 
 ! 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              
+!           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
       END DO   
+
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU   ,  'U', -1., zspgV   , 'V', -1. )
+      CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_ai,  'U', -1., ztauy_ai, 'V', -1. )
 
       CALL iom_put( 'zmsk00x'    , zmsk00x  )   ! MV DEBUG
@@ -431,6 +436,5 @@
          ! In the outer loop, one needs to update all RHS terms
          ! with explicit velocity dependencies (viscosities, coriolis, ocean stress)
-         ! as a function of uc
-         !
+         ! as a function of "current" velocities (uc, vc)
       
          !------------------------------------------
@@ -451,5 +455,6 @@
 
          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 ( i_out == nn_nout_vp ) CALL iom_put( 'zds'          , zds          )   ! MV DEBUG
 
          IF( lwp )   WRITE(numout,*) ' outer loop  1a i_out : ', i_out
@@ -464,7 +469,7 @@
 
                ! shear**2 at T points (doc eq. A16)
-               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
-                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
-                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+               zds2  = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &    + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &    ) * 0.25_wp * r1_e1e2t(ji,jj)
               
                ! divergence at T points
@@ -475,14 +480,14 @@
                
                ! tension at T points
-               zdt  = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-                  &   - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
-                  &   ) * r1_e1e2t(ji,jj)
-               zdt2 = zdt * zdt
+               zdt   = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &    - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &    ) * r1_e1e2t(ji,jj)
+               zdt2  = zdt * zdt
                
                ! delta at T points
-               zdeltat = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+               zdeltat(ji,jj)        = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
                
                ! delta* at T points (following Lemieux and Dupont, GMD 2020)
-               zdeltastar_t(ji,jj) = zdeltat + rn_creepl
+               zdeltastar_t(ji,jj)   = zdeltat(ji,jj) + rn_creepl
               
                ! P/delta at T-points
@@ -490,15 +495,19 @@
                
                ! Temporary zzt and zet factors at T-points
-               zzt(ji,jj)     = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)
-               zet(ji,jj)     = zzt(ji,jj)     * z1_ecc2 
+               zzt(ji,jj)            = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)
+               zet(ji,jj)            = zzt(ji,jj)     * z1_ecc2 
                           
             END DO
          END DO
-         
-         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
+
+!        zzt(:,:) = 0. ! test
+!        zet(:,:) = 0. ! test
+         
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. , zdeltat, 'T', 1.)
+
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zet'        , zet      )   ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zdeltat'       , zdeltat        ) ! MV DEBUG
 
          IF( lwp )   WRITE(numout,*) ' outer loop  1b i_out : ', i_out
@@ -515,7 +524,9 @@
             END DO
          END DO
+
+!        zef(:,:) = 0. ! test
          
          CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. )
-         CALL iom_put( 'zef'           , zef            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zef'           , zef            ) ! MV DEBUG
          IF( lwp )   WRITE(numout,*) ' outer loop  1c i_out : ', i_out
 
@@ -551,5 +562,5 @@
                            &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_c(ji,jj  ) + e2u(ji-1,jj  ) * zu_c(ji-1,jj  ) )  &
                            &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) )
-         
+
              END DO
          END DO
@@ -560,8 +571,8 @@
          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
-         CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG
-         CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG
 
          IF( lwp )   WRITE(numout,*) ' outer loop  1f i_out : ', i_out
@@ -575,9 +586,10 @@
 
          ! --- Stress contributions at T-points         
-         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
+         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1 & zs2
             DO ji = 2, jpi ! 
             
                ! 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) - 1.0_wp )
+               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)
                                             
                ! sig2 contribution to RHS of U-equation at T-points            
@@ -585,5 +597,6 @@
 
                ! 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) - 1.0_wp )
+               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)
 
                ! sig2 contribution to RHS of V-equation  at T-points
@@ -593,12 +606,15 @@
          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
+         ! TEST
+         ! zs2_rhsv(:,:) = 0._wp
+
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG
          
          ! a priori, no lbclnk, because rhsu is only used in the inner domain
          
-         ! --- Stress contributions at f-points         
+         ! --- Stress contributions at F-points         
          ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing
          ! My guess is that this is the way to enforce boundary conditions on strain rate tensor
@@ -610,17 +626,20 @@
                
                ! sig12 contribution to RHS of U equation at F-points 
-               zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) - r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj)
+               zs12_rhsu(ji,jj) =   zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) + r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj)
                
                ! sig12 contribution to RHS of V equation at F-points
-               zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) - r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj)
+               zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) + r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj)
 
             END DO
          END DO
+         !
+         ! TEST
+         ! zs12_rhsv(:,:) = 0._wp
 
          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
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG
 
          ! a priori, no lbclnk, because rhsu are only used in the inner domain
@@ -631,4 +650,5 @@
          DO jj = 2, jpj - 1
             DO ji = 2, jpi - 1               
+
                ! --- U component of internal force contribution to RHS at U points
                zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & 
@@ -640,12 +660,12 @@
                zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * &
                   &           (    e1v(ji,jj)    * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) )                                                                 &
-                  &      +         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj)     * zs2_rhsv(ji,jj) )     &
+                  &      -         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj)   * e1t(ji,jj)   * zs2_rhsv(ji,jj) )     &
                   &      + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj)   * e2f(ji,jj)   * zs12_rhsv(ji,jj)  - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) )
-                  
+
             END DO
          END DO
 
-         CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG
-         CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG
          
          !---------------------------
@@ -657,7 +677,6 @@
             DO ji = 2, jpi - 1
             
-               ! still miss ice ocean stress and acceleration contribution
                zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj)
-               zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsu(ji,jj)
+               zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsv(ji,jj)
 
             END DO
@@ -668,10 +687,10 @@
          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
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG
          
          ! inner domain calculations -> no lbclnk
@@ -712,22 +731,22 @@
                zfac3      = 2._wp * zfac * r1_e1u(ji,jj)
 
-               zt12U      = - zfac1 * zzt(ji+1,jj)
                zt11U      =   zfac1 * zzt(ji,jj)
-         
-               zt22U      = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj)
+               zt12U      =   zfac1 * zzt(ji+1,jj)
+         
                zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)
-         
-               zt122U     = - zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)
+               zt22U      =   zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj)
+         
                zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)
+               zt122U     =   zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)
                
                !
                ! Linear system coefficients
                !
-               zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj)
-               zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj)
-               zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj)
-         
-               zDU(ji,jj) =   zt121U * r1_e1u(ji,jj-1)
-               zEU(ji,jj) =   zt122U * r1_e1u(ji,jj+1)
+               zAU(ji,jj) = -   zt11U           * e2u(ji-1,jj) -   zt21U           * r1_e2u(ji-1,jj)
+               zBU(ji,jj) =   ( zt11U + zt12U ) * e2u(ji,jj)   + ( zt21U + zt22U ) * r1_e2u(ji,jj)   + ( zt121U + zt122U ) * r1_e1u(ji,jj)
+               zCU(ji,jj) = -   zt12U           * e2u(ji+1,jj) -   zt22U           * r1_e2u(ji+1,jj)
+         
+               zDU(ji,jj) =     zt121U * r1_e1u(ji,jj-1)
+               zEU(ji,jj) =     zt122U * r1_e1u(ji,jj+1)
               
                !
@@ -737,26 +756,26 @@
                !
                zfac       = 0.5_wp * r1_e1e2v(ji,jj)
-               zfac1      =         zfac * e2v(ji,jj)
+               zfac1      =         zfac * e1v(ji,jj)
                zfac2      =         zfac * r1_e1v(ji,jj)
                zfac3      = 2._wp * zfac * r1_e2v(ji,jj)
          
-               zt12V      = - zfac1 * zzt(ji,jj+1)
                zt11V      =   zfac1 * zzt(ji,jj)
-         
+               zt12V      =   zfac1 * zzt(ji,jj+1)
+         
+               zt21V      =   zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)
                zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1)
-               zt21V      = - zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)
-         
+         
+               zt121V     =   zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)
                zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)
-               zt121V     = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)
-         
+
                !
                ! Linear system coefficients
                !
-               zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1)
-               zBV(ji,jj) =  ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj)
-               zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1)
-         
-               zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here
-               zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj)
+               zAV(ji,jj) = -   zt11V           * e1v(ji,jj-1) -   zt21V           * r1_e1v(ji,jj-1)
+               zBV(ji,jj) =   ( zt11V + zt12V ) * e1v(ji,jj)   + ( zt21V + zt22V ) * r1_e1v(ji,jj)   + ( zt122V + zt121V ) * r1_e2v(ji,jj)
+               zCV(ji,jj) = -   zt12V           * e1v(ji,jj+1) -   zt22V           * r1_e1v(ji,jj+1)
+         
+               zDV(ji,jj) =     zt121V * r1_e2v(ji-1,jj)
+               zEV(ji,jj) =     zt122V * r1_e2v(ji+1,jj)
                   
                !-----------------------------------------------------
@@ -764,10 +783,10 @@
                !-----------------------------------------------------
                zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj)
-               zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj)
+               zBV(ji,jj) = zBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj)
          
             END DO
          END DO
 
-         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU  , 'U', 1., zAV  , 'V',  1. )
+         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. )
@@ -775,14 +794,15 @@
          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
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG
+
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG
+         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG
 
       !------------------------------------------------------------------------------!
@@ -808,6 +828,6 @@
             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
+            !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
@@ -864,4 +884,8 @@
                      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)
+
                         DO ji = 3, jpi - 1
 
@@ -962,5 +986,9 @@
                      DO ji = ji_min, jpi - 1, nn_zebra_vp 
                      
-                        DO jj = 3, jpj - 1
+                        ! 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)
+
+                        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 ) )
@@ -1137,21 +1165,23 @@
 
       ! 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. )
-
-      IF ( lwp ) WRITE(numout,*) ' Velocity replaced '
+      ! info: probably wrong because when doing du/dt = div ( sigma ) we have spurious waves instead of u = v = 0
+
+!     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. )
+
+!     IF ( lwp ) WRITE(numout,*) ' Velocity replaced '
 
       ! END DEBUG
