Index: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90	(revision 13564)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90	(revision 13566)
@@ -172,6 +172,6 @@
          END DO
       END DO
-      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. )
-      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1. )
       !
       !
@@ -392,4 +392,17 @@
          ENDIF
          !
+         ! --- Lateral boundary conditions --- !
+         IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN
+            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
+               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp )
+         ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN
+            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
+               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp )
+         ELSE
+            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp )
+         ENDIF
+         CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp )
+         CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp )
+         !
          !== Open water area ==!
          zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
@@ -400,6 +413,5 @@
             END DO
          END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. )
-         !
+         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp )
          !
          ! --- Ensure non-negative fields and in-bound thicknesses --- !
@@ -461,6 +473,5 @@
       !!             work on H (and not V). It is partly related to the multi-category approach
       !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
-      !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter
-      !!             since sv_i and e_i are still good.
+      !!             concentration is small). We also limit S and T.
       !!----------------------------------------------------------------------
       REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0)
@@ -506,5 +517,5 @@
       IF( pamsk == 0._wp ) THEN
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 2, jpjm1
                DO ji = 1, fs_jpim1
                   IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
@@ -516,4 +527,8 @@
                   ENDIF
                   !
+               END DO
+            END DO
+            DO jj = 1, jpjm1
+               DO ji = fs_2, fs_jpim1
                   IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
                      zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
@@ -553,8 +568,14 @@
       IF( PRESENT( pua_ho ) ) THEN
          DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = 1, fs_jpim1
+                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl)
+                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl)
+              END DO
+            END DO
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
-                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
+               DO ji = fs_2, fs_jpim1
+                  pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
+                  pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
               END DO
             END DO
@@ -573,5 +594,4 @@
          END DO
       END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. )
       !
    END SUBROUTINE adv_umx
@@ -614,5 +634,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
+               DO jj = 1, jpj
                   DO ji = 1, fs_jpim1
                      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)
@@ -622,5 +642,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO jj = 2, jpjm1
+               DO jj = 1, jpj
                   DO ji = fs_2, fs_jpim1
                      ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              &
@@ -631,9 +651,8 @@
                END DO
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
             !
             DO jl = 1, jpl              !-- flux in y-direction
                DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
+                  DO ji = fs_2, fs_jpim1
                      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 DO
@@ -645,5 +664,5 @@
             DO jl = 1, jpl              !-- flux in y-direction
                DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
+                  DO ji = 1, jpi
                      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 DO
@@ -653,5 +672,5 @@
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
                DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
+                  DO ji = 1, jpi
                      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)
@@ -661,8 +680,7 @@
                END DO
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
+               DO jj = 2, jpjm1
                   DO ji = 1, fs_jpim1
                      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)
@@ -717,7 +735,11 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 1, jpj
                DO ji = 1, fs_jpim1
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
+               END DO
+            END DO
+            DO jj = 1, jpjm1
+               DO ji = 1, jpi
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
                END DO
@@ -737,5 +759,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
+               DO jj = 1, jpj
                   DO ji = 1, fs_jpim1
                      pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
@@ -746,5 +768,5 @@
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO jj = 2, jpjm1
+               DO jj = 1, jpj
                   DO ji = fs_2, fs_jpim1
                      ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              &
@@ -755,9 +777,8 @@
                END DO
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
 
             DO jl = 1, jpl              !-- flux in y-direction
                DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
+                  DO ji = fs_2, fs_jpim1
                      pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
                   END DO
@@ -770,5 +791,5 @@
             DO jl = 1, jpl              !-- flux in y-direction
                DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
+                  DO ji = 1, jpi
                      pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
                   END DO
@@ -779,5 +800,5 @@
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
                DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
+                  DO ji = 1, jpi
                      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)
@@ -787,8 +808,7 @@
                END DO
             END DO
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO jj = 1, jpjm1
+               DO jj = 2, jpjm1
                   DO ji = 1, fs_jpim1
                      pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
@@ -953,5 +973,5 @@
          !        
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 2, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
@@ -964,5 +984,5 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 2, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -976,5 +996,5 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 2, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -992,5 +1012,5 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 2, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -1008,5 +1028,5 @@
          !
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 2, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -1031,5 +1051,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO jj = 1, jpjm1
+            DO jj = 2, jpjm1
                DO ji = 1, fs_jpim1
                   IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
@@ -1043,5 +1063,5 @@
       !                                                     !-- High order flux in i-direction  --!
       DO jl = 1, jpl
-         DO jj = 1, jpjm1
+         DO jj = 2, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
@@ -1112,5 +1132,5 @@
          DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
+               DO ji = fs_2, fs_jpim1
                   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) ) )
@@ -1122,5 +1142,5 @@
          DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
+               DO ji = fs_2, fs_jpim1
                   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)   &
@@ -1133,5 +1153,5 @@
          DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
+               DO ji = fs_2, fs_jpim1
                   zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                   zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1148,5 +1168,5 @@
          DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
+               DO ji = fs_2, fs_jpim1
                   zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                   zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1163,5 +1183,5 @@
          DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
+               DO ji = fs_2, fs_jpim1
                   zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                   zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1186,5 +1206,5 @@
          DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
+               DO ji = fs_2, fs_jpim1
                   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) )  &
@@ -1198,5 +1218,5 @@
       DO jl = 1, jpl
          DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
+            DO ji = fs_2, fs_jpim1
                pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
             END DO
@@ -1235,7 +1255,11 @@
       ! --------------------------------------------------
       DO jl = 1, jpl
-         DO jj = 1, jpjm1
+         DO jj = 2, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
+            END DO
+         END DO
+         DO jj = 1, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
             END DO
@@ -1352,5 +1376,5 @@
       ! ---------------------------------
       DO jl = 1, jpl
-         DO jj = 1, jpjm1
+         DO jj = 2, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
@@ -1366,5 +1390,5 @@
 
          DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
                zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
