Index: NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90	(revision 10925)
+++ NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90	(revision 10929)
@@ -36,29 +36,25 @@
    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
    
-   ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
-   INTEGER ::   kn_limiter = 1
-
-   ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),
-   !    interpolated T at u/v points can be non-zero while it should
-   !    (because of the high order of the advection scheme). Thus set it to 0 in this case
-   LOGICAL ::   ll_icedge = .TRUE.
-
+   ! advection for S and T:    dVS/dt = -div( uA * uHS / u ) => ll_advS = F
+   !                        or dVS/dt = -div( uV * uS  / u ) => ll_advS = T
+   LOGICAL ::   ll_advS = .FALSE.
+   !
+   ! alternate directions for upstream
+   LOGICAL ::   ll_upsxy = .TRUE.
+   !
+   ! alternate directions for high order
+   LOGICAL ::   ll_hoxy = .TRUE.
+   !
    ! if T interpolated at u/v points is negative or v_i < 1.e-6
    !    then interpolate T at u/v points using the upstream scheme
    LOGICAL ::   ll_neg = .TRUE.
-   
-   ! alternate directions for upstream
-   LOGICAL ::   ll_upsxy = .TRUE.
-
-   ! alternate directions for high order
-   LOGICAL ::   ll_hoxy = .TRUE.
-   
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zmski_small, zmskj_small
+   !
+   ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
+   INTEGER ::   kn_limiter = 1
+   !
    ! prelimiter: use it to avoid overshoot in H
    LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
-
-   ! advection for S and T:    dVS/dt = -div( uA * uHS / u ) => ll_advS = F
-   !                        or dVS/dt = -div( uV * uS  / u ) => ll_advS = T
-   LOGICAL ::   ll_advS = .FALSE.
-
+   !
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -102,5 +98,5 @@
       INTEGER  ::   icycle                  ! number of sub-timestep for the advection
       REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
-      REAL(wp) ::   zdt
+      REAL(wp) ::   zdt, zvi_cen
       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication
       REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
@@ -134,5 +130,5 @@
          END DO
       END DO
-      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
+      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
       !
       !
@@ -185,4 +181,26 @@
          ELSEWHERE                        ;   z1_aip(:,:,:) = 0.
          END WHERE
+         !
+         ! setup a mask where advection will be upstream
+         IF( ll_neg ) THEN
+            IF( .NOT. ALLOCATED(zmski_small) )   ALLOCATE( zmski_small(jpi,jpj,jpl) ) 
+            IF( .NOT. ALLOCATED(zmskj_small) )   ALLOCATE( zmskj_small(jpi,jpj,jpl) ) 
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
+                     IF( zvi_cen < epsi06) THEN ; zmski_small(ji,jj,jl) = 0.
+                     ELSE                       ; zmski_small(ji,jj,jl) = 1. ; ENDIF
+                     zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
+                     IF( zvi_cen < epsi06) THEN ; zmskj_small(ji,jj,jl) = 0.
+                     ELSE                       ; zmskj_small(ji,jj,jl) = 1. ; ENDIF
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         !
+         ! ----------------------- !
+         ! ==> start advection <== !
+         ! ----------------------- !
          !
          ! set u*A=u for advection of A only 
@@ -319,5 +337,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. )
          !
          !
@@ -326,9 +344,9 @@
          ! Remove negative values (conservation is ensured)
          !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
-         CALL ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
          !
          ! Make sure ice thickness is not too big
          !    (because ice thickness can be too large where ice concentration is very small)
-         CALL Hbig( zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
 
       END DO
@@ -762,7 +780,4 @@
                      &                                                                                        * pamsk           &
                      &                             ) * pdt ) * tmask(ji,jj,1)
-                  !!clem test
-                  !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) )
-                  !!clem test
                END DO
             END DO
@@ -794,7 +809,4 @@
                      &                                                                                        * pamsk           &
                      &                             ) * pdt ) * tmask(ji,jj,1) 
-                  !!clem test
-                  !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) )
-                  !!clem test
                END DO
             END DO
@@ -837,5 +849,5 @@
       !
       INTEGER  ::   ji, jj, jl             ! dummy loop indices
-      REAL(wp) ::   zcu, zdx2, zdx4, zvi_cen2        !   -      -
+      REAL(wp) ::   zcu, zdx2, zdx4        !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4
       !!----------------------------------------------------------------------
@@ -948,19 +960,4 @@
       END SELECT
       !
-      ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),
-      !    interpolated T at u/v points can be non-zero while it should
-      !    (because of the high order of the advection scheme). Thus set it to 0 in this case
-      IF( ll_icedge ) THEN
-         DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  IF( pt(ji,jj,jl) <= 0._wp .AND. pu(ji,jj) >= 0._wp ) THEN
-                     pt_u(ji,jj,jl) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-      !
       ! if pt at u-point is negative then use the upstream value
       !    this should not be necessary if a proper sea-ice mask is set in Ultimate
@@ -970,6 +967,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zvi_cen2 = 0.5_wp * ( v_i(ji+1,jj,jl) + v_i(ji,jj,jl) )
-                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN
+                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zmski_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)   &
                         &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -1010,5 +1006,5 @@
       !
       INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zcv, zdy2, zdy4, zvi_cen2    !   -      -
+      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4
       !!----------------------------------------------------------------------
@@ -1118,19 +1114,4 @@
       END SELECT
       !
-      ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),
-      !    interpolated T at u/v points can be non-zero while it should
-      !    (because of the high order of the advection scheme). Thus set it to 0 in this case
-      IF( ll_icedge ) THEN
-         DO jl = 1, jpl
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  IF( pt(ji,jj,jl) <= 0._wp .AND. pv(ji,jj) >= 0._wp ) THEN
-                     pt_v(ji,jj,jl) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-      !
       ! if pt at v-point is negative then use the upstream value
       !    this should not be necessary if a proper sea-ice mask is set in Ultimate
@@ -1140,6 +1121,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zvi_cen2 = 0.5_wp * ( v_i(ji,jj+1,jl) + v_i(ji,jj,jl) )
-                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN
+                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zmskj_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) )  &
                         &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1528,5 +1508,5 @@
 
 
-   SUBROUTINE Hbig( phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE Hbig  ***
@@ -1545,4 +1525,5 @@
       !! ** input   : Max thickness of the surrounding 9-points
       !!-------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step
       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip
@@ -1551,8 +1532,9 @@
       !
       INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
-      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra
+      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra
       REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
       !!-------------------------------------------------------------------
       !
+      z1_dt = 1._wp / pdt
       !
       DO jl = 1, jpl
@@ -1584,6 +1566,6 @@
                      zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
                      !
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
                      !
                      pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
@@ -1598,6 +1580,6 @@
                   IF( zvs_excess > 0._wp ) THEN
                      zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
                      !
                      pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
Index: NEMO/releases/release-4.0/src/ICE/icevar.F90
===================================================================
--- NEMO/releases/release-4.0/src/ICE/icevar.F90	(revision 10925)
+++ NEMO/releases/release-4.0/src/ICE/icevar.F90	(revision 10929)
@@ -229,5 +229,5 @@
                   IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
                      !
-                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
+                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i             ! Energy of melting e(S,T) [J.m-3]
                      ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C]
                      ! Conversion q(S,T) -> T (second order equation)
@@ -236,5 +236,5 @@
                      t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
                      !
-                  ELSE                                !--- no ice
+                  ELSE                                   !--- no ice
                      t_i(ji,jj,jk,jl) = rt0
                   ENDIF
@@ -537,5 +537,5 @@
 
 
-   SUBROUTINE ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_zapneg ***
@@ -543,6 +543,5 @@
       !! ** Purpose :   Remove negative sea ice fields and correct fluxes
       !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
-      !
+      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
       REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
@@ -555,6 +554,10 @@
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
-      !!-------------------------------------------------------------------
-      !
+      !
+      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
+      REAL(wp) ::   z1_dt
+      !!-------------------------------------------------------------------
+      !
+      z1_dt = 1._wp / pdt
       !
       DO jl = 1, jpl       !==  loop over the categories  ==!
@@ -567,5 +570,5 @@
                DO ji = 1 , jpi
                   IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 >0
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
                      pe_i(ji,jj,jk,jl) = 0._wp
                   ENDIF
@@ -578,5 +581,5 @@
                DO ji = 1 , jpi
                   IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
                      pe_s(ji,jj,jk,jl) = 0._wp
                   ENDIF
@@ -591,13 +594,13 @@
             DO ji = 1 , jpi
                IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
                   pv_i   (ji,jj,jl) = 0._wp
                ENDIF
                IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
                   pv_s   (ji,jj,jl) = 0._wp
                ENDIF
                IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
-                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice
+                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
                   psv_i  (ji,jj,jl) = 0._wp
                ENDIF
