Index: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5034)
+++ branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5600)
@@ -102,5 +102,5 @@
       !!                 and charge ellipse.
       !!                 The user should make sure that the parameters
-      !!                 nevp, telast and creepl maintain stress state
+      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
       !!                 on the charge ellipse for plastic flow
       !!                 e.g. in the Canadian Archipelago
@@ -108,5 +108,4 @@
       !! References : Hunke and Dukowicz, JPO97
       !!              Bouillon et al., Ocean Modelling 2009
-      !!              Vancoppenolle et al., Ocean Modelling 2008
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
@@ -117,18 +116,17 @@
       CHARACTER (len=50) ::   charout
       REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
-      REAL(wp) ::   za, zstms, zmask   ! local scalars
-      REAL(wp) ::   zc1, zc2, zc3             ! ice mass
-
-      REAL(wp) ::   dtevp              ! time step for subcycling
-      REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity
-      REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
-      REAL(wp) ::   zu_ice2, zv_ice1   !
-      REAL(wp) ::   zddc, zdtc         ! delta on corners and on centre
-      REAL(wp) ::   zdst               ! shear at the center of the grid point
-      REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface
-      REAL(wp) ::   sigma1, sigma2     ! internal ice stress
+      REAL(wp) ::   za, zstms          ! local scalars
+      REAL(wp) ::   zc1, zc2, zc3      ! ice mass
+
+      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling
+      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity
+      REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars
+      REAL(wp) ::   zu_ice2, zv_ice1              !
+      REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre
+      REAL(wp) ::   zdst                          ! shear at the center of the grid point
+      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface
+      REAL(wp) ::   sigma1, sigma2                ! internal ice stress
 
       REAL(wp) ::   zresm         ! Maximal error on ice velocity
-      REAL(wp) ::   zdummy        ! dummy argument
       REAL(wp) ::   zintb, zintn  ! dummy argument
 
@@ -139,8 +137,9 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                           
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                           
+      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points
       REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points
       
       REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
@@ -152,8 +151,11 @@
                                                               !   ocean surface (ssh_m) if ice is not embedded
                                                               !   ice top surface if ice is embedded   
+
+      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter
+      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
       !!-------------------------------------------------------------------
 
       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
+      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
       CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
       CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
@@ -161,14 +163,14 @@
 #if  defined key_lim2 && ! defined key_lim2_vp
 # if defined key_agrif
-     USE ice_2, vt_s => hsnm
-     USE ice_2, vt_i => hicm
+      USE ice_2, vt_s => hsnm
+      USE ice_2, vt_i => hicm
 # else
-     vt_s => hsnm
-     vt_i => hicm
+      vt_s => hsnm
+      vt_i => hicm
 # endif
-     at_i(:,:) = 1. - frld(:,:)
+      at_i(:,:) = 1. - frld(:,:)
 #endif
 #if defined key_agrif && defined key_lim2 
-    CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
+      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
 #endif
       !
@@ -186,19 +188,17 @@
 
 #if defined key_lim3
-      CALL lim_itd_me_icestrength( ridge_scheme_swi )      ! LIM-3: Ice strength on T-points
-#endif
-
-!CDIR NOVERRCHK
+      CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points
+#endif
+
       DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
-!CDIR NOVERRCHK
          DO ji = 1 , jpi
 #if defined key_lim3
-            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj)
+            zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj)
 #endif
 #if defined key_lim2
-            zpresh(ji,jj) = tms(ji,jj) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
-#endif
-            ! tmi = 1 where there is ice or on land
-            tmi(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj)
+            zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
+#endif
+            ! zmask = 1 where there is ice or on land
+            zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)
          END DO
       END DO
@@ -206,20 +206,13 @@
       ! Ice strength on grid cell corners (zpreshc)
       ! needed for calculation of shear stress 
-!CDIR NOVERRCHK
       DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
          DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
-            zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
-               &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
-               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &
-               &              tms(ji,jj)     * wght(ji+1,jj+1,1,1)
-            zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
-               &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
-               &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
-               &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   &
-               &             ) / MAX( zstms, epsd )
+            zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   &
+               &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1)
+            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
+               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     &
+               &             ) / MAX( zstms, zepsi )
          END DO
       END DO
-
       CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
       !
@@ -236,21 +229,19 @@
       !  zcorl2: Coriolis parameter on V-points                            
       !  (ztagnx,ztagny): wind stress on U/V points                       
-      !  u_oce1: ocean u component on u points                           
       !  v_oce1: ocean v component on u points                          
       !  u_oce2: ocean u component on v points                         
-      !  v_oce2: ocean v component on v points                        
 
       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
-          !                                            
-          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
-          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
+         !                                            
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
          zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
-          !
-          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
-          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         !
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
-          !
+         !
          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
-          !
+         !
       ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
          zpice(:,:) = ssh_m(:,:)
@@ -260,36 +251,33 @@
          DO ji = fs_2, fs_jpim1
 
-            zc1 = tms(ji  ,jj  ) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
-            zc2 = tms(ji+1,jj  ) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
-            zc3 = tms(ji  ,jj+1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
-
-            zt11 = tms(ji  ,jj) * e1t(ji  ,jj)
-            zt12 = tms(ji+1,jj) * e1t(ji+1,jj)
-            zt21 = tms(ji,jj  ) * e2t(ji,jj  )
-            zt22 = tms(ji,jj+1) * e2t(ji,jj+1)
+            zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
+            zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
+            zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
+
+            zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj)
+            zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj)
+            zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  )
+            zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1)
 
             ! Leads area.
-            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
-            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
+            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi )
+            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi )
 
             ! Mass, coriolis coeff. and currents
-            zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+epsd)
-            zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+epsd)
-            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   &
-               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd )
-            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) )   &
-               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
+            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
+            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
+            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   &
+               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi )
+            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   &
+               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi )
             !
-            u_oce1(ji,jj)  = u_oce(ji,jj)
-            v_oce2(ji,jj)  = v_oce(ji,jj)
-
             ! Ocean has no slip boundary condition
-            v_oce1(ji,jj)  = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)    &
-               &                 +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) &
-               &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)  
-
-            u_oce2(ji,jj)  = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)     &
-               &                 +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) &
-               &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
+            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      &
+               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  &
+               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
+
+            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      &
+               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  &
+               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
 
             ! Wind stress at U,V-point
@@ -303,6 +291,6 @@
             ! include it later
 
-            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
-            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
+            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
 
             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
@@ -318,9 +306,14 @@
       !
       ! Time step for subcycling
-      dtevp  = rdt_ice / nevp
+      dtevp  = rdt_ice / nn_nevp
+#if defined key_lim3
+      dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )
+#else
       dtotel = dtevp / ( 2._wp * telast )
-
+#endif
+      z1_dtotel = 1._wp / ( 1._wp + dtotel )
+      z1_dtevp  = 1._wp / dtevp
       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
-      ecc2 = ecc * ecc
+      ecc2 = rn_ecc * rn_ecc
       ecci = 1. / ecc2
 
@@ -331,5 +324,5 @@
 
       !                                               !----------------------!
-      DO jter = 1 , nevp                              !    loop over jter    !
+      DO jter = 1 , nn_nevp                           !    loop over jter    !
          !                                            !----------------------!        
          DO jj = k_j1, k_jpj-1
@@ -339,5 +332,5 @@
 
          DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
+            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask
 
                !  
@@ -360,122 +353,80 @@
                !
                !
-               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
-                  &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
-                  &             +e1v(ji,jj)*v_ice(ji,jj)                      &
-                  &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
-                  &             )                                             &
-                  &            / area(ji,jj)
-
-               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
-                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
-                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
-                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
-                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
-                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
-                  &         )                                              &
-                  &        / area(ji,jj)
+               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &            ) * r1_e12t(ji,jj)
+
+               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &         ) * r1_e12t(ji,jj)
 
                !
-               zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                &
-                  &            -u_ice(ji,jj)/e1u(ji,jj)                    &
-                  &           )*e1f(ji,jj)*e1f(ji,jj)                      &
-                  &          +( v_ice(ji+1,jj)/e2v(ji+1,jj)                &
-                  &            -v_ice(ji,jj)/e2v(ji,jj)                    &
-                  &           )*e2f(ji,jj)*e2f(ji,jj)                      &
-                  &         )                                              &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
-                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
-                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
-
-
-               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
-                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
-                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
-
-               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
-                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
-                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
-
-            END DO
-         END DO
-         CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond.
-
-!CDIR NOVERRCHK
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+                  &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   &
+                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
+
+
+               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
+                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
+
+               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
+                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond.
+         
          DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1
 
                !- Calculate Delta at centre of grid cells
-               zdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
-                  &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          &
-                  &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          &
-                  &          - e1v(ji, jj-1) * u_ice2(ji,jj-1)          &
-                  &          )                                          &
-                  &         / area(ji,jj)
-
-               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
-               delta_i(ji,jj) = delta + creepl
-               !-Calculate stress tensor components zs1 and zs2 
-               !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) )  &
-                  &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
-               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) )  &
-                  &         / ( 1._wp + dtotel )
-
-            END DO
-         END DO
-
-         CALL lbc_lnk( zs1(:,:), 'T', 1. )
-         CALL lbc_lnk( zs2(:,:), 'T', 1. )
-
-!CDIR NOVERRCHK
-         DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
-            DO ji = fs_2, fs_jpim1
+               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
+                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
+                  &            ) * r1_e12t(ji,jj)
+
+               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
+               delta_i(ji,jj) = delta + rn_creepl
+
                !- Calculate Delta on corners
-               zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
-                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
-                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
-                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
-                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
-                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
-                  &         )                                               &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
-
-               zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
-                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
-                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
-                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
-                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
-                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
-                  &         )                                               &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
-
-               zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
-
-               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
-                  &          ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) )  &
-                  &          / ( 1.0 + dtotel ) 
-
-            END DO ! ji
-         END DO ! jj
-
-         CALL lbc_lnk( zs12(:,:), 'F', 1. )
-
+               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &    ) * r1_e12f(ji,jj)
+
+               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &    ) * r1_e12f(ji,jj)
+
+               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl
+
+               !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide).
+               zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   &
+                  &          ) * z1_dtotel
+               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   &
+                  &          ) * z1_dtotel
+               !-Calculate stress tensor component zs12 at corners
+               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  &
+                  &          ) * z1_dtotel 
+
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+ 
          ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
          DO jj = k_j1+1, k_jpj-1
             DO ji = fs_2, fs_jpim1
                !- contribution of zs1, zs2 and zs12 to zf1
-               zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) &
-                  &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) &
-                  &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) &
-                  &             ) / ( e1u(ji,jj)*e2u(ji,jj) )
+               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
+                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  &
+                  &                ) * r1_e12u(ji,jj)
                ! contribution of zs1, zs2 and zs12 to zf2
-               zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) &
-                  &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) &
-                  &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    &
-                  zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) &
-                  &             ) / ( e1v(ji,jj)*e2v(ji,jj) )
+               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
+                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  &
+                  &               )  * r1_e12v(ji,jj)
             END DO
          END DO
@@ -487,23 +438,19 @@
          IF (MOD(jter,2).eq.0) THEN 
 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
-                  z0           = zmass1(ji,jj)/dtevp
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
+                  z0           = zmass1(ji,jj) * z1_dtevp
 
                   ! SB modif because ocean has no slip boundary condition
-                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         &
-                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
-                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
-                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
-                     (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
-                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
-                     za*(u_oce1(ji,jj))
-                  zcca         = z0+za
+                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     &
+                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
+                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
+                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
+                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
+                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
-
+                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
                END DO
             END DO
@@ -511,5 +458,5 @@
             CALL lbc_lnk( u_ice(:,:), 'U', -1. )
 #if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nevp, 'U' )
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
 #endif
 #if defined key_bdy
@@ -517,23 +464,19 @@
 #endif         
 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
 
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
-                  z0           = zmass2(ji,jj)/dtevp
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
+                  z0           = zmass2(ji,jj) * z1_dtevp
                   ! SB modif because ocean has no slip boundary condition
-                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     &
-                     &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
-                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
-                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
-                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
-                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
-                     za2ct(ji,jj) + za*(v_oce2(ji,jj))
-                  zcca         = z0+za
+                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
+                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
+                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
+                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
+                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) )
+                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
-
+                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
                END DO
             END DO
@@ -541,5 +484,5 @@
             CALL lbc_lnk( v_ice(:,:), 'V', -1. )
 #if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nevp, 'V' )
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
 #endif
 #if defined key_bdy
@@ -548,23 +491,19 @@
 
          ELSE 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
-                  z0           = zmass2(ji,jj)/dtevp
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
+                  z0           = zmass2(ji,jj) * z1_dtevp
                   ! SB modif because ocean has no slip boundary condition
-                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      &
-                     &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
-                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)   
-
-                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &
-                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
-                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
-                     za2ct(ji,jj) + za*(v_oce2(ji,jj))
-                  zcca         = z0+za
+                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
+                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
+                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)   
+
+                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &
+                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )
+                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
-
+                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
                END DO
             END DO
@@ -572,5 +511,5 @@
             CALL lbc_lnk( v_ice(:,:), 'V', -1. )
 #if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nevp, 'V' )
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
 #endif
 #if defined key_bdy
@@ -578,27 +517,24 @@
 #endif         
 
-!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
                DO ji = fs_2, fs_jpim1
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
-                  z0           = zmass1(ji,jj)/dtevp
-                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      &
-                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
-                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
-
-                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
-                     (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
-                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
-                     za*(u_oce1(ji,jj))
-                  zcca         = z0+za
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
+                  z0           = zmass1(ji,jj) * z1_dtevp
+                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       &
+                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
+                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
+
+                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
+                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
+                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
-               END DO ! ji
-            END DO ! jj
+                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
+               END DO
+            END DO
 
             CALL lbc_lnk( u_ice(:,:), 'U', -1. )
 #if defined key_agrif && defined key_lim2
-            CALL agrif_rhg_lim2( jter, nevp, 'U' )
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
 #endif
 #if defined key_bdy
@@ -611,8 +547,7 @@
             !---  Convergence test.
             DO jj = k_j1+1 , k_jpj-1
-               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
-                  ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
-            END DO
-            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
+               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
+            END DO
+            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
             IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
          ENDIF
@@ -625,24 +560,20 @@
       ! 4) Prevent ice velocities when the ice is thin
       !------------------------------------------------------------------------------!
-      ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the
-      ! ocean velocity, 
-      ! This prevents high velocity when ice is thin
-!CDIR NOVERRCHK
+      ! If the ice volume is below zvmin then ice velocity should equal the
+      ! ocean velocity. This prevents high velocity when ice is thin
       DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1
-            zdummy = vt_i(ji,jj)
-            IF ( zdummy .LE. hminrhg ) THEN
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
                u_ice(ji,jj) = u_oce(ji,jj)
                v_ice(ji,jj) = v_oce(ji,jj)
-            ENDIF ! zdummy
+            ENDIF
          END DO
       END DO
 
-      CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
-      CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
+      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. )
+
 #if defined key_agrif && defined key_lim2
-      CALL agrif_rhg_lim2( nevp , nevp, 'U' )
-      CALL agrif_rhg_lim2( nevp , nevp, 'V' )
+      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' )
+      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' )
 #endif
 #if defined key_bdy
@@ -653,93 +584,64 @@
       DO jj = k_j1+1, k_jpj-1 
          DO ji = fs_2, fs_jpim1
-            zdummy = vt_i(ji,jj)
-            IF ( zdummy .LE. hminrhg ) THEN
-               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
-                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
-                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
-
-               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
-                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
-                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
-            ENDIF ! zdummy
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
+               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     &
+                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
+                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
+
+               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
+                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
+            ENDIF 
          END DO
       END DO
 
-      CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 
-      CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
+      CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. )
 
       ! Recompute delta, shear and div, inputs for mechanical redistribution 
-!CDIR NOVERRCHK
       DO jj = k_j1+1, k_jpj-1
-!CDIR NOVERRCHK
-         DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
+         DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask
             !- divu_i(:,:), zdt(:,:): divergence and tension at centre 
             !- zds(:,:): shear on northeast corner of grid cells
-            zdummy = vt_i(ji,jj)
-            IF ( zdummy .LE. hminrhg ) THEN
-
-               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
-                  &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
-                  &             +e1v(ji,jj)*v_ice(ji,jj)                      &
-                  &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
-                  &            )                                              &
-                  &            / area(ji,jj)
-
-               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
-                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
-                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
-                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
-                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
-                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
-                  &         )                                              &
-                  &        / area(ji,jj)
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
+
+               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   &
+                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   &
+                  &            ) * r1_e12t(ji,jj)
+
+               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
+                  &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
+                  &         ) * r1_e12t(ji,jj)
                !
                ! SB modif because ocean has no slip boundary condition 
-               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1)              &
-                  &           - u_ice(ji,jj)   / e1u(ji,jj) )              &
-                  &           * e1f(ji,jj) * e1f(ji,jj)                    &
-                  &          + ( v_ice(ji+1,jj) / e2v(ji+1,jj)             &
-                  &            - v_ice(ji,jj)  / e2v(ji,jj) )              &
-                  &           * e2f(ji,jj) * e2f(ji,jj) )                  &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
-                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
-                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
-
-               zdst = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    &
-                  &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    &
-                  &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    &
-                  &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj)
-
-               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
-               delta_i(ji,jj) = delta + creepl
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     &
+                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
+
+               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    &
+                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj)
+
+               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
+               delta_i(ji,jj) = delta + rn_creepl
             
-            ENDIF ! zdummy
-
-         END DO !jj
-      END DO !ji
+            ENDIF
+         END DO
+      END DO
       !
       !------------------------------------------------------------------------------!
       ! 5) Store stress tensor and its invariants
       !------------------------------------------------------------------------------!
-      !
       ! * Invariants of the stress tensor are required for limitd_me
       !   (accelerates convergence and improves stability)
       DO jj = k_j1+1, k_jpj-1
          DO ji = fs_2, fs_jpim1
-            ! begin TECLIM change 
-            zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
-               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &   
-               &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &   
-               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 
+            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &   
+               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 
             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
-            ! end TECLIM change
          END DO
       END DO
 
       ! Lateral boundary condition
-      CALL lbc_lnk( divu_i (:,:), 'T', 1. )
-      CALL lbc_lnk( delta_i(:,:), 'T', 1. )
-      ! CALL lbc_lnk( shear_i(:,:), 'F', 1. )
-      CALL lbc_lnk( shear_i(:,:), 'T', 1. )
+      CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. )
 
       ! * Store the stress tensor for the next time step
@@ -772,5 +674,5 @@
             DO jj = k_j1+1, k_jpj-1
                DO ji = 2, jpim1
-                  IF (zpresh(ji,jj) .GT. 1.0) THEN
+                  IF (zpresh(ji,jj) > 1.0) THEN
                      sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
                      sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
@@ -786,5 +688,5 @@
       !
       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
+      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
       CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
       CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
