Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90	(revision 13536)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90	(revision 13544)
@@ -17,4 +17,5 @@
    USE ice            ! sea-ice: variables
    USE icedyn_rhg_evp ! sea-ice: EVP rheology
+   USE icedyn_rhg_vp  ! sea-ice: VP  rheology
    USE icectl         ! sea-ice: control prints
    !
@@ -82,7 +83,4 @@
          !                             !---------------------------------!
          CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
-         !         
-      END SELECT
-
       !                                !---------------------------------------------!
       CASE( np_rhgVP  )                ! Viscous-Plastic rheology                    !
@@ -166,5 +164,4 @@
       !
       IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
-      IF( ln_rhg_VP   )   CALL rhg_lsr_rst( 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE ice_dyn_rhg_init
Index: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90	(revision 13536)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90	(revision 13544)
@@ -169,5 +169,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
       !! --- diags
-      REAL(wp) ::   zsig1, zsig2, zsig12
+      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i         
       !! --- SIMIP diags
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 13536)
+++ NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90	(revision 13544)
@@ -1,33 +1,2 @@
-! TODOLOIST
-!
-! Commit
-! Compile
-!
-! - Re-calculate internal force diagnostic (currently incorrect, see end of the routine)
-! - QC:compare code to mitGCM
-! - QC: comparing EVP and VP codes 
-!
-! Clarify
-! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ?
-! --- Swap of independent term in the U-V solvers ? -> JF Lemieux says maybe
-! --- Is stress tensor for restart needed ? No!!!
-!
-! Test
-! --- Try ADI for u-v solver
-!
-! Add
-! - Charge ellipse as an output (good one from Lemieux and Dupont)
-! - Bottom drag (landfast ice)
-! - Tensile strength
-! - agrif
-! - bdy
-!
-! Write documentation
-!
-! Optimize
-! - subroutinize common parts (diagnostics)
-! - namelist: rationalize common parameters EVP/VP
-
-
 MODULE icedyn_rhg_vp
    !!======================================================================
@@ -73,5 +42,5 @@
 
    !! for convergence tests
-   INTEGER ::   ncvgid   ! netcdf file id
+   INTEGER ::   ncvgid        ! netcdf file id
    INTEGER ::   nvarid_ucvg   ! netcdf variable id
    INTEGER ::   nvarid_ures
@@ -82,4 +51,5 @@
    INTEGER ::   nvarid_mke
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
+
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -157,5 +127,5 @@
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
       !!
-      LOGICAL ::   ll_u_iterate, ll_viterate   ! continue iteration or not
+      LOGICAL ::   ll_u_iterate, ll_v_iterate   ! continue iteration or not
       !
       INTEGER ::   ji, jj, jn          ! dummy loop indices
@@ -171,5 +141,6 @@
       REAL(wp) ::   zm2, zm3, zmassU, zmassV                            ! ice/snow mass and volume
       REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars
-      REAL(wp) ::   zp_deltastar_f
+      REAL(wp) ::   zp_deltastar_f                                      ! 
+      REAL(wp) ::   zu_cV, zv_cU                                        ! 
       REAL(wp) ::   zfac, zfac1, zfac2, zfac3
       REAL(wp) ::   zt12U, zt11U, zt22U, zt21U, zt122U, zt121U
@@ -178,5 +149,5 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
-      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   za_iU  , za_iV                      ! ice fraction on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! Acceleration term contribution to RHS
       REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step
@@ -191,5 +162,4 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_c, zv_c                      ! "current" ice velocity (m/s), average of previous two OL iterates
-      REAL(wp), DIMENSION(jpi,jpj) ::   zu_cV, zv_cU                    ! "current" u (v) ice velocity at V (U) point (m/s)
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_b, zv_b                      ! velocity at previous sub-iterate
       REAL(wp), DIMENSION(jpi,jpj) ::   zuerr, zverr                    ! absolute U/Vvelocity difference between current and previous sub-iterates
@@ -209,5 +179,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zrhsu, zrhsv                    ! U and V RHS 
       REAL(wp), DIMENSION(jpi,jpj) ::   zAU, zBU, zCU, zDU, zEU         ! Linear system coefficients, U equation
-      REAL(wp), DIMENSION(jpi,jpj) ::   zAV, zBV, zCV, zDV, zEV, zFV    ! Linear system coefficients, V equation
+      REAL(wp), DIMENSION(jpi,jpj) ::   zAV, zBV, zCV, zDV, zEV         ! Linear system coefficients, V equation
       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
@@ -222,13 +192,11 @@
       REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
       !! --- diags
-      REAL(wp) ::   zsig1, zsig2, zsig12
+      REAL(wp) ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f, zten_i   ! stress tensor components
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s, SIMIP)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s, SIMIP)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s, SIMIP)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s, SIMIP)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s, SIMIP)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s, SIMIP)      
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice, zdiag_ymtrp_ice ! X/Y-component of ice mass transport (kg/s, SIMIP)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw, zdiag_ymtrp_snw ! X/Y-component of snow mass transport (kg/s, SIMIP)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp, zdiag_yatrp         ! X/Y-component of area transport (m2/s, SIMIP)
                         
       !!----------------------------------------------------------------------------------------------------------------------
@@ -251,7 +219,7 @@
       END DO
       
-      IF ( ln_zebra_vp ) THEN; nn_nzebra_vp = 2; ELSE; nn_nzebra_vp = 1; ENDIF 
-      
-      nn_nvp = nn_out_vp * nn_inn_vp ! maximum number of iterations
+      IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2; ELSE; nn_zebra_vp = 1; ENDIF 
+      
+      nn_nvp = nn_nout_vp * nn_ninn_vp ! maximum number of iterations
       
       zrhoco = rau0 * rn_cio 
@@ -341,6 +309,6 @@
 
             ! Ice fraction at U-V points
-            zaU(ji,jj)      = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
-            zaV(ji,jj)      = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+            za_iU(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            za_iV(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
 
             ! Snow and ice mass at U-V points
@@ -368,6 +336,6 @@
             
             ! Wind stress
-            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
-            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+            ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj)
+            ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj)
 
             ! Force due to sea surface tilt(- m*g*GRAD(ssh))
@@ -380,7 +348,7 @@
 
             ! Mask for lots of ice (1) or little ice (0)
-            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+            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. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
+            IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
             ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
 
@@ -499,8 +467,8 @@
                 
                 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
-                zCwU(ji,jj)          = zaU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  &
-                  &                                              + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
-                zCwV(ji,jj)          = zaV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  &
-                  &                                              + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
+                zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  &
+                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
+                zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  &
+                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
                  
                 !--- Ocean-ice drag contributions to RHS 
@@ -577,5 +545,5 @@
                               (    e2u(ji,jj)    * ( zs1_rhsu(ji+1,jj) - zs1_rhsu(ji,jj) )                                                                 &
                   &      +         r1_e2u(ji,jj) * ( e2t(ji+1,jj) * e2t(ji+1,jj) * zs2_rhsu(ji+1,jj) - e2t(ji,jj)   * e2t(ji,jj)   * zs2_rhsu(ji,jj) )     &
-                  &      + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj)   * e1f(ji,jj)   * zs12_rhsu(ji,jj)  - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) ) 
+                  &      + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj)   * e1f(ji,jj)   * zs12_rhsu(ji,jj)  - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) ) )
                   
                ! --- V component of internal force contribution to RHS at V points
@@ -583,5 +551,5 @@
                   &           (    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) )     &
-                  &      + 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) )
+                  &      + 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
@@ -723,5 +691,5 @@
                   ! what follows could be subroutinized...
                   
-                  DO jn = 1, nn_nzebra_vp ! "zebra" loop (! red-black-sor!!! )
+                  DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! )
                   
                      ! OPT: could be even better optimized with a true red-black SOR
@@ -735,5 +703,5 @@
                      zBU_prime(:,:) = 0._wp
          
-                     DO jj = jj_min, jpj - 1, nn_nzebra_vp
+                     DO jj = jj_min, jpj - 1, nn_zebra_vp
 
                         !-------------------------------------------
@@ -756,8 +724,8 @@
      
                            ! right hand side
-                           zFU(ji,jj) = ( zrhsu(ji,jj) &                                 ! right-hand side terms
-                               &      +   zAA3                                           ! boundary condition translation
-                               &      +   zDU(ji,jj) * u_ice(ji,jj-1)                    ! internal force, j-1
-                               &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1
+                           zFU(ji,jj) = ( zrhsu(ji,jj) &                                   ! right-hand side terms
+                               &      +   zAA3                                           & ! boundary condition translation
+                               &      +   zDU(ji,jj) * u_ice(ji,jj-1)                    & ! internal force, j-1
+                               &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1)   ! internal force, j+1
                     
                         END DO
@@ -806,5 +774,5 @@
                   ! MV OPT: what follows could be subroutinized...
                   
-                  DO jn = 1, nn_nzebra_vp ! "zebra" loop
+                  DO jn = 1, nn_zebra_vp ! "zebra" loop
       
                      IF ( jn == 1 ) THEN   ;   ji_min = 2 
@@ -816,5 +784,5 @@
                      zBV_prime(:,:) = 0._wp
 
-                     DO ji = ji_min, jpi - 1, nn_nzebra_vp 
+                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
                         
                         !!! It is intentional to have a ji then jj loop for V-velocity
@@ -836,8 +804,8 @@
      
                            ! right hand side
-                           zFV(ji,jj) = ( zrhsv(ji,jj) &              ! right-hand side terms
-                               &        + zAA3                        ! boundary condition translation
-                               &        + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1
-                               &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1
+                           zFV(ji,jj) = ( zrhsv(ji,jj) &                                   ! right-hand side terms
+                               &        + zAA3                                           & ! boundary condition translation
+                               &        + zDV(ji,jj) * v_ice(ji-1,jj)                    & ! internal force, j-1
+                               &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1)   ! internal force, j+1
                     
                         END DO
@@ -855,5 +823,5 @@
                         v_ice(ji,jpj-1)  = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1)  ! last row
                         DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed
-                           v_ice(ji,jj)   = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj) / zBV_prime(ji,jj)
+                           v_ice(ji,jj)   = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) / zBV_prime(ji,jj)
                         END DO            
             
@@ -902,6 +870,6 @@
                   
                   ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
-                  IF ( i_inn > 1 && zuerr_max > rn_uerr_max_vp ) THEN
-                      IF ( lwp ) " VP rheology error was too large : ", zu_err_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
+                  IF ( i_inn > 1 .AND. zuerr_max > rn_uerr_max_vp ) THEN
+                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
                       ll_u_iterate = .FALSE.
                   ENDIF
@@ -909,5 +877,5 @@
                   ! - Stop if error small enough
                   IF ( zuerr_max < rn_uerr_min_vp ) THEN                                        
-                      IF ( lwp ) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! "
+                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! "
                       ll_u_iterate = .FALSE.
                   ENDIF
@@ -933,12 +901,12 @@
                   
                   ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
-                  IF ( i_inn > 1 && zverr_max > rn_uerr_max_vp ) THEN
-                      IF ( lwp ) " VP rheology error was too large : ", zv_err_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
+                  IF ( i_inn > 1 .AND. zverr_max > rn_uerr_max_vp ) THEN
+                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zverr_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
                       ll_v_iterate = .FALSE.
                   ENDIF
                   
                   ! - Stop if error small enough
-                  IF ( zverr_max < rn_verr_min ) THEN                                        
-                      IF ( lwp ) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! "
+                  IF ( zverr_max < rn_uerr_min_vp ) THEN                                        
+                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! "
                       ll_v_iterate = .FALSE.
                   ENDIF
@@ -965,5 +933,5 @@
       !------------------------------------------------------------------------------!
 
-      END ! End of outer loop (i_out) =============================================================================================
+      END DO ! End of outer loop (i_out) =============================================================================================
 
       !------------------------------------------------------------------------------!
@@ -1077,7 +1045,31 @@
       IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
          & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
-         !
-         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
-            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+
+         ALLOCATE( ztaux_oi(jpi,jpj) , ztauy_oi(jpi,jpj) )
+
+         !--- Recalculate oceanic stress at last inner iteration
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
+
+                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
+                zu_cV            = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
+                zv_cU            = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
+                
+                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
+                zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( u_ice(ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) )  &
+                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
+                zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( v_ice(ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) )  &
+                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
+                 
+                !--- Ocean-ice stress
+                ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                
+            END DO
+         END DO
+         
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
+!            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
          !
          CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
@@ -1085,6 +1077,8 @@
          CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
          CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
-         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
-         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
+!        CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
+!        CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
+
+         DEALLOCATE( ztaux_oi , ztauy_oi )
 
       ENDIF
@@ -1167,18 +1161,31 @@
       IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
          & iom_use('corstrx') .OR. iom_use('corstry') ) THEN
+
+         ! --- Recalculate Coriolis stress at last inner iteration
+         DO jj = 2, jpj - 1
+            DO ji = 2, jpi - 1
+                ! --- U-component 
+                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  &
+                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  &
+                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+            END DO
+         END DO
          !
          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
             &                                 zCorU, 'U', -1., zCorV, 'V', -1. )
-
+         !
          CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
          CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
          CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
          CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+
       ENDIF
       
       IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
 
-
-         ! Recalculate internal forces (divergence of stress tensor)
+         ! Recalculate internal forces (divergence of stress tensor) at last inner iteration
          DO jj = 2, jpj - 1
             DO ji = 2, jpi - 1
@@ -1291,9 +1298,8 @@
       INTEGER           ::   it, idtime, istatus
       INTEGER           ::   ji, jj          ! dummy loop indices
-      REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zmke, zu, zv ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)   ::   zu_res(:,:), zv_res(:,:), zvel2(:,:) ! local arrays
+      REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays
                                                                              
       CHARACTER(len=20) ::   clname
-      ::   zres           ! check convergence
       !!----------------------------------------------------------------------
 
@@ -1352,16 +1358,16 @@
       DO jj = 2, jpj - 1
          DO ji = 2, jpi - 1                                      
-            zu_res(ji,jj)  = zrhsu(ji,jj) + zDU(ji,jj) * pu(ji,jj-1) + zEU(ji,jj) * pu(ji,jj+1) &
-               &           - zAU(ji,jj) * pu(ji-1,jj) - zBU(ji,jj) * pu(ji,jj) - zCU(ji,jj) * pu(ji+1,jj)
+            zu_res(ji,jj)  = prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) &
+               &           - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj)
                            
-            zv_res(ji,jj)  = zrhsv(ji,jj) + zDV(ji,jj) * pv(ji-1,jj) + zEV(ji,jj) * pv(ji+1,jj) &
-               &           - zAV(ji,jj) * pv(ji,jj-1) - zBV(ji,jj) * pv(ji,jj) - zCV(ji,jj) * pv(ji,jj+1)                           
+            zv_res(ji,jj)  = prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) &
+               &           - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1)                           
           END DO
        END DO                  
-       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:) )
-       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:) )
-       zu_res_mean = SQRT( zu_resmean / pglob_area )
-       zv_res_mean = SQRT( zv_resmean / pglob_area )
-       zvelres     = MEAN( zu_res_mean, zv_res_mean )
+       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:,1) )
+       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:,1) )
+       zu_res_mean = SQRT( zu_res_mean / pglob_area )
+       zv_res_mean = SQRT( zv_res_mean / pglob_area )
+       zvelres     = 0.5_wp * ( zu_res_mean + zv_res_mean )
                                           
        ! -- Global mean kinetic energy per unit area (J/m2)  
@@ -1374,5 +1380,5 @@
        END DO
        
-       zmke = 0.5_wp * globsum( 'ice_rhg_vp', zmt(:,:) * e1e2t(:,:) * zvel2(:,:) )
+       zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) )
                   
          !                                                ! ==================== !
