Index: /NEMO/trunk/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 15048)
+++ /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 15049)
@@ -268,10 +268,10 @@
             &                          , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  &
             &                          , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume
-            &                          , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  )
-         CALL lbc_lnk( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity
+            &                          , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  &
+            &                          , z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity
             &                          , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  &
             &                          , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration
-            &                          , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  )
-         CALL lbc_lnk( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age
+            &                          , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  &
+            &                          , z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age
             &                          , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  )
          CALL lbc_lnk( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy
@@ -280,11 +280,16 @@
             &                          , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  )
          IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
-            CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
-               &                          , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
-               &                          , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
-               &                          , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )
-            IF ( ln_pnd_lids ) THEN
-               CALL lbc_lnk( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume
-                  &                          , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  )
+            IF( ln_pnd_lids ) THEN
+               CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
+                  &                          , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
+                  &                          , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
+                  &                          , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  &
+                  &                          , z0vl , 'T', 1._wp, sxvl , 'T', -1._wp, syvl , 'T', -1._wp  & ! melt pond lid volume
+                  &                          , sxxvl, 'T', 1._wp, syyvl, 'T',  1._wp, sxyvl, 'T',  1._wp  )
+            ELSE
+               CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
+                  &                          , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
+                  &                          , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
+                  &                          , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )
             ENDIF
          ENDIF
@@ -766,5 +771,5 @@
       !
       DO jl = 1, jpl
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
@@ -813,5 +818,5 @@
       !                                           ! -- check e_i/v_i -- !
       DO jl = 1, jpl
-         DO_3D( 1, 1, 1, 1, 1, nlay_i )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -827,5 +832,5 @@
       !                                           ! -- check e_s/v_s -- !
       DO jl = 1, jpl
-         DO_3D( 1, 1, 1, 1, 1, nlay_s )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s )
             IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
                ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -870,5 +875,5 @@
       ! -- check snow load -- !
       DO jl = 1, jpl
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
@@ -1196,7 +1201,6 @@
          END_2D
       END DO
-      
    END SUBROUTINE icemax3D
-         
+
    SUBROUTINE icemax4D( pice , pmax )
       !!---------------------------------------------------------------------
@@ -1237,5 +1241,4 @@
          END DO
       END DO
-
    END SUBROUTINE icemax4D
 
Index: /NEMO/trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 15048)
+++ /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 15049)
@@ -164,10 +164,11 @@
       !
       ! --- define velocity for advection: u*grad(H) --- !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( nn_hls-1, nn_hls, nn_hls, nn_hls )
          IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
          ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
          ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
          ENDIF
-
+      END_2D
+      DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls )
          IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
          ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
@@ -204,8 +205,10 @@
             IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) )
             DO jl = 1, jpl
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 0, nn_hls, nn_hls )
                   zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
                   IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0
                   ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF
+               END_2D
+               DO_2D( nn_hls, nn_hls, 1, 0 )
                   zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
                   IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0
@@ -583,5 +586,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
                pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
                pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
@@ -594,5 +597,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 1, 1 )
+               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )
                   pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
                END_2D
@@ -600,5 +603,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D( 0, 0, 1, 1 )
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls )
                   ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              &
                      &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -609,5 +612,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 0, 0, 1, 0 )
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
                   pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl)
                END_2D
@@ -617,5 +620,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 1, 1, 0 )
+               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )
                   pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
                END_2D
@@ -623,5 +626,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D( 1, 1, 0, 0 )
+               DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 )
                   ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  &
                      &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -632,5 +635,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 0, 0 )
+               DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
                   pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl)
                END_2D
@@ -642,5 +645,5 @@
       !
       DO jl = 1, jpl                    !-- after tracer with upstream scheme
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   &
                &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) &
@@ -651,5 +654,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1.0_wp )
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1.0_wp )
 
    END SUBROUTINE upstream
@@ -681,8 +684,8 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 1 )
+            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )
                pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
             END_2D
-            DO_2D( 1, 1, 1, 0 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )
                pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
             END_2D
@@ -701,5 +704,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 1, 1 )
+               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
                END_2D
@@ -708,5 +711,5 @@
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D( 0, 0, 1, 1 )
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls )
                   ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              &
                      &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -717,5 +720,5 @@
 
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 0, 0, 1, 0 )
+               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
                END_2D
@@ -726,5 +729,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 1, 1, 0 )
+               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
                END_2D
@@ -733,5 +736,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D( 1, 1, 0, 0 )
+               DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 )
                   ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  &
                      &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -742,5 +745,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 0, 0 )
+               DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
                END_2D
@@ -785,10 +788,10 @@
          !
          !                                                        !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
+         CALL ultimate_x( nn_hls, pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
          !                                                        !--  limiter in x --!
          IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, nn_hls, nn_hls )
                zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pubox(ji,jj   ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t  (ji,jj) &
                   &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) &
@@ -797,11 +800,10 @@
             END_2D
          END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
          !
          !                                                        !--  ultimate interpolation of pt at v-point  --!
          IF( ll_hoxy ) THEN
-            CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
+            CALL ultimate_y( 0, pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
          ELSE
-            CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
+            CALL ultimate_y( 0, pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
          ENDIF
          !                                                        !--  limiter in y --!
@@ -812,10 +814,10 @@
          !
          !                                                        !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
+         CALL ultimate_y( nn_hls, pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
          !                                                        !--  limiter in y --!
          IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
          !                                                        !--  advective form update in zpt  --!
          DO jl = 1, jpl
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls, nn_hls, 0, 0 )
                zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pvbox(ji,jj   ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t  (ji,jj) &
                   &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) &
@@ -824,11 +826,10 @@
             END_2D
          END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
          !
          !                                                        !--  ultimate interpolation of pt at u-point  --!
          IF( ll_hoxy ) THEN
-            CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
+            CALL ultimate_x( 0, pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
          ELSE
-            CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
+            CALL ultimate_x( 0, pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
          ENDIF
          !                                                        !--  limiter in x --!
@@ -842,5 +843,5 @@
 
 
-   SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
+   SUBROUTINE ultimate_x( kloop, pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -852,4 +853,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
       !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kloop     ! either 0 or nn_hls depending on the order of the call
       REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0)
       INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
@@ -867,23 +869,40 @@
       !                                                     !--  Laplacian in i-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 0, 0 )      ! First derivative (gradient)
+         DO_2D( nn_hls, nn_hls-1, kloop, kloop )      ! First derivative (gradient)
             ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
          END_2D
-         DO_2D( 0, 0, 0, 0 )      ! Second derivative (Laplacian)
+         DO_2D( nn_hls-1, nn_hls-1, kloop, kloop )    ! Second derivative (Laplacian)
             ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
          END_2D
-      END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1.0_wp )
+!!$         DO jj = 2, jpjm1         ! First derivative (gradient)
+!!$            DO ji = 1, jpim1
+!!$               ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+!!$            END DO
+!!$            !                     ! Second derivative (Laplacian)
+!!$            DO ji = 2, jpim1
+!!$               ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
+!!$            END DO
+!!$         END DO
+      END DO
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1.0_wp )
       !
       !                                                     !--  BiLaplacian in i-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 0, 0 )      ! Third derivative
+         DO_2D( 1, 0, kloop, kloop )                  ! Third derivative
             ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
          END_2D
-         DO_2D( 0, 0, 0, 0 )      ! Fourth derivative
-               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
-         END_2D
-      END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1.0_wp )
+         DO_2D( 0, 0, kloop, kloop )                  ! Fourth derivative
+            ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
+         END_2D
+!!$         DO jj = 2, jpjm1         ! Third derivative
+!!$            DO ji = 1, jpim1
+!!$               ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+!!$            END DO
+!!$            !                     ! Fourth derivative
+!!$            DO ji = 2, jpim1
+!!$               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
+!!$            END DO
+!!$         END DO
+      END DO
       !
       !
@@ -893,5 +912,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 1, 0, kloop, kloop )
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
                   &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -902,5 +921,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 1, 0, kloop, kloop )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
@@ -912,5 +931,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 1, 0, kloop, kloop )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -926,5 +945,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 1, 0, kloop, kloop )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -939,6 +958,8 @@
       CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
          !
-         DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+         CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1.0_wp )
+         !
+         DO jl = 1, jpl
+            DO_2D( 1, 0, kloop, kloop )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -961,5 +982,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D( 1, 0, 0, 0 )
+            DO_2D( 1, 0, kloop, kloop )
                IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
@@ -979,5 +1000,5 @@
 
 
-   SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
+   SUBROUTINE ultimate_y( kloop, pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_y  ***
@@ -989,4 +1010,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
       !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kloop     ! either 0 or nn_hls depending on the order of the call
       REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0)
       INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
@@ -1004,23 +1026,22 @@
       !                                                     !--  Laplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 0, 0, 1, 0 )         ! First derivative (gradient)
+         DO_2D( kloop, kloop, nn_hls, nn_hls-1 )      ! First derivative (gradient)
             ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
          END_2D
-         DO_2D( 0, 0, 0, 0 )         ! Second derivative (Laplacian)
+         DO_2D( kloop, kloop, nn_hls-1, nn_hls-1 )    ! Second derivative (Laplacian)
             ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1.0_wp )
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1.0_wp )
       !
       !                                                     !--  BiLaplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 0, 0, 1, 0 )         ! Third derivative
+         DO_2D( kloop, kloop, 1, 0 )                  ! Third derivative
             ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
          END_2D
-         DO_2D( 0, 0, 0, 0 )         ! Fourth derivative
+         DO_2D( kloop, kloop, 0, 0 )                  ! Fourth derivative
             ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1.0_wp )
       !
       !
@@ -1029,5 +1050,5 @@
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( kloop, kloop, 1, 0 )
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
                   &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1037,5 +1058,5 @@
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( kloop, kloop, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
@@ -1046,5 +1067,5 @@
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( kloop, kloop, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1059,5 +1080,5 @@
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( kloop, kloop, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1071,6 +1092,9 @@
          !
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
-         DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+         !
+         CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1.0_wp )
+         !
+         DO jl = 1, jpl
+            DO_2D( kloop, kloop, 1, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1093,5 +1117,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D( 0, 0, 1, 0 )
+            DO_2D( kloop, kloop, 1, 0 )
                IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
@@ -1293,12 +1317,12 @@
       !
       DO jl = 1, jpl
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls, nn_hls-1, 0, 0 )
             zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.0_wp)   ! lateral boundary cond.
-
-      DO jl = 1, jpl
-         DO_2D( 0, 0, 0, 0 )
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.0_wp)   ! lateral boundary cond.
+
+      DO jl = 1, jpl
+         DO_2D( nn_hls-1, 0, 0, 0 )
             uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
 
@@ -1361,5 +1385,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp)   ! lateral boundary cond.
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp)   ! lateral boundary cond.
       !
    END SUBROUTINE limiter_x
@@ -1384,12 +1408,12 @@
       !
       DO jl = 1, jpl
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, nn_hls, nn_hls-1 )
             zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.0_wp)   ! lateral boundary cond.
-
-      DO jl = 1, jpl
-         DO_2D( 0, 0, 0, 0 )
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.0_wp)   ! lateral boundary cond.
+
+      DO jl = 1, jpl
+         DO_2D( 0, 0, nn_hls-1, 0 )
             vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
 
@@ -1453,5 +1477,5 @@
          END_2D
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.0_wp)   ! lateral boundary cond.
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.0_wp)   ! lateral boundary cond.
       !
    END SUBROUTINE limiter_y
@@ -1488,5 +1512,5 @@
       !
       DO jl = 1, jpl
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
@@ -1535,5 +1559,5 @@
       !                                           ! -- check e_i/v_i -- !
       DO jl = 1, jpl
-         DO_3D( 1, 1, 1, 1, 1, nlay_i )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -1549,5 +1573,5 @@
       !                                           ! -- check e_s/v_s -- !
       DO jl = 1, jpl
-         DO_3D( 1, 1, 1, 1, 1, nlay_s )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s )
             IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
                ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
@@ -1592,5 +1616,5 @@
       ! -- check snow load -- !
       DO jl = 1, jpl
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
@@ -1649,5 +1673,4 @@
          END_2D
       END DO
-
    END SUBROUTINE icemax3D
 
@@ -1690,5 +1713,4 @@
          END DO
       END DO
-      
    END SUBROUTINE icemax4D
 
Index: /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 15048)
+++ /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 15049)
@@ -184,9 +184,9 @@
       !
       ! for diagnostics and convergence tests
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
       END_2D
       IF( nn_rhg_chkcvg > 0 ) THEN
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
          END_2D
@@ -259,5 +259,11 @@
       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
 
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zm1          = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  ! Ice/snow mass at U-V points
+         zmf  (ji,jj) = zm1 * ff_t(ji,jj)                            ! Coriolis at T points (m*f)
+         zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin )                   ! dt/m at T points (for alpha and beta coefficients)
+      END_2D
+      
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
 
          ! ice fraction at U-V points
@@ -276,10 +282,4 @@
          u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
 
-         ! Coriolis at T points (m*f)
-         zmf(ji,jj)      = zm1 * ff_t(ji,jj)
-
-         ! dt/m at T points (for alpha and beta coefficients)
-         zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
-
          ! m/dt
          zmU_t(ji,jj)    = zmassU * z1_dtevp
@@ -305,10 +305,9 @@
 
       END_2D
-      CALL lbc_lnk( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
       !
       !                                  !== Landfast ice parameterization ==!
       !
       IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ! ice thickness at U-V points
             zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
@@ -327,5 +326,5 @@
          !
       ELSE                               !-- no landfast
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ztaux_base(ji,jj) = 0._wp
             ztauy_base(ji,jj) = 0._wp
@@ -351,5 +350,5 @@
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
 
             ! shear at F points
@@ -386,9 +385,9 @@
 
          ! P/delta at T points
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
          END_2D
 
-         DO_2D( 0, 1, 0, 1 )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
+         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
 
             ! divergence at T points (duplication to avoid communications)
@@ -425,10 +424,10 @@
          ! Save beta at T-points for further computations
          IF( ln_aEVP ) THEN
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
             END_2D
          ENDIF
 
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
 
             ! alpha for aEVP
@@ -450,5 +449,5 @@
 
          ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             !                   !--- U points
             zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
@@ -478,5 +477,5 @@
          IF( MOD(jter,2) == 0 ) THEN ! even iterations
             !
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                !                 !--- tau_io/(v_oce - v_ice)
                zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
@@ -522,11 +521,5 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
-            !
-#if defined key_agrif
-!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
-            CALL agrif_interp_ice( 'V' )
-#endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+            IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
             !
             DO_2D( 0, 0, 0, 0 )
@@ -574,15 +567,11 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
-            !
-#if defined key_agrif
-!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
-            CALL agrif_interp_ice( 'U' )
-#endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            IF( nn_hls == 1 ) THEN   ;   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
+            ELSE                     ;   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
+            ENDIF
             !
          ELSE ! odd iterations
             !
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                !                 !--- tau_io/(u_oce - u_ice)
                zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
@@ -628,11 +617,5 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
-            !
-#if defined key_agrif
-!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
-            CALL agrif_interp_ice( 'U' )
-#endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
             !
             DO_2D( 0, 0, 0, 0 )
@@ -680,14 +663,19 @@
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
-            !
-#if defined key_agrif
-!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
-            CALL agrif_interp_ice( 'V' )
-#endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+            IF( nn_hls == 1 ) THEN   ;   CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
+            ELSE                     ;   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
+            ENDIF
             !
          ENDIF
-
+         !
+#if defined key_agrif
+!!       CALL agrif_interp_ice( 'U', jter, nn_nevp )
+!!       CALL agrif_interp_ice( 'V', jter, nn_nevp )
+         CALL agrif_interp_ice( 'U' )
+         CALL agrif_interp_ice( 'V' )
+#endif
+         IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+         IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+         !
          ! convergence test
          IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
@@ -701,5 +689,5 @@
       ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
       !------------------------------------------------------------------------------!
-      DO_2D( 1, 0, 1, 0 )
+      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
 
          ! shear at F points
@@ -777,5 +765,5 @@
          ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
             ! Ice stresses
@@ -810,5 +798,5 @@
          ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
             ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
