Index: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icedyn_adv_umx.F90	(revision 10744)
+++ NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icedyn_adv_umx.F90	(revision 12065)
@@ -11,8 +11,8 @@
    !!   'key_si3'                                       SI3 sea-ice model
    !!----------------------------------------------------------------------
-   !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme
+   !!   ice_dyn_adv_umx   : update the tracer fields
    !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
-   !!   macho             : ???
-   !!   nonosc_ice        : compute monotonic tracer fluxes by a non-oscillatory algorithm 
+   !!   macho             : compute the fluxes
+   !!   nonosc_ice        : limit the fluxes using a non-oscillatory algorithm 
    !!----------------------------------------------------------------------
    USE phycst         ! physical constant
@@ -23,4 +23,5 @@
    !
    USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
    USE lib_mpp        ! MPP library
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
@@ -31,27 +32,22 @@
 
    PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
-      
-   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
-   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
-   
-   ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
-   INTEGER ::   kn_limiter = 1
-
-   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
-   !   clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why)
-   !         but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration 
-   LOGICAL ::   ll_neg = .TRUE.
-   
-   ! alternate directions for upstream
-   LOGICAL ::   ll_upsxy = .TRUE.
-
-   ! alternate directions for high order
-   LOGICAL ::   ll_hoxy = .TRUE.
-   
-   ! prelimiter: use it to avoid overshoot in H
-   !   clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why)
-   LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
-
-
+   !
+   INTEGER, PARAMETER ::   np_advS = 1         ! advection for S and T:    dVS/dt = -div(      uVS     ) => np_advS = 1
+   !                                                                    or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
+   !                                                                    or dVS/dt = -div( uV * uS  / u ) => np_advS = 3
+   INTEGER, PARAMETER ::   np_limiter = 1      ! limiter: 1 = nonosc
+   !                                                      2 = superbee
+   !                                                      3 = h3
+   LOGICAL            ::   ll_upsxy  = .TRUE.   ! alternate directions for upstream
+   LOGICAL            ::   ll_hoxy   = .TRUE.   ! alternate directions for high order
+   LOGICAL            ::   ll_neg    = .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_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
+   !
+   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6
+   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120
+   !
+   INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small
+   !
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -63,5 +59,5 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
+   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
@@ -78,4 +74,7 @@
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness
       REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
@@ -84,5 +83,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
@@ -92,18 +91,43 @@
       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), DIMENSION(1)           ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
+      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
       REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
       !
-      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
-      !     When needed, the advection split is applied at the next time-step in order to avoid blocking global comm.
-      !     ...this should not affect too much the stability... Was ok on the tests we did...
+      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
+                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
+                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
+                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
+               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
+                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
+                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
+                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
+               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
+                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
+                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
+                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
+      !
+      !
+      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
+      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
+      !              this should not affect too much the stability
       zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
       zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
@@ -115,5 +139,4 @@
       ELSE                         ;   icycle = 1
       ENDIF
-      
       zdt = rdt_ice / REAL(icycle)
 
@@ -121,5 +144,11 @@
       zudy(:,:) = pu_ice(:,:) * e2u(:,:)
       zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
-
+      !
+      ! setup transport for each ice cat 
+      DO jl = 1, jpl
+         zu_cat(:,:,jl) = zudy(:,:)
+         zv_cat(:,:,jl) = zvdx(:,:)
+      END DO
+      !
       ! --- define velocity for advection: u*grad(H) --- !
       DO jj = 2, jpjm1
@@ -153,60 +182,178 @@
          END WHERE
          !
-         ! set u*a=u for advection of A only 
-         DO jl = 1, jpl
-            zua_ho(:,:,jl) = zudy(:,:)
-            zva_ho(:,:,jl) = zvdx(:,:)
-         END DO
-         
+         ! setup a mask where advection will be upstream
+         IF( ll_neg ) THEN
+            IF( .NOT. ALLOCATED(imsk_small) )   ALLOCATE( imsk_small(jpi,jpj,jpl) ) 
+            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_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   ;   imsk_small(ji,jj,jl) = 0
+                     ELSE                         ;   imsk_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   ;   jmsk_small(ji,jj,jl) = 0
+                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         !
+         ! ----------------------- !
+         ! ==> start advection <== !
+         ! ----------------------- !
+         !
+         !== Ice area ==!
          zamsk = 1._wp
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area
-         zamsk = 0._wp
-         !
-         zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                ! Ice volume
-         !
-         zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                ! Snw volume
-         !
-         zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               ! Salt content
-         !
-         zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )               ! Age content
-         !
-         DO jk = 1, nlay_i
-            zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) )   ! Ice heat content
-         END DO
-         !
-         DO jk = 1, nlay_s
-            zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) )   ! Snw heat content
-         END DO
-         !
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, &
+            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
+         !
+         !                             ! --------------------------------- !
+         IF( np_advS == 1 ) THEN       ! -- advection form: -div( uVS ) -- !
+            !                          ! --------------------------------- !
+            zamsk = 0._wp
+            !== Ice volume ==!
+            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_i, zua_ups, zva_ups )
+            !== Snw volume ==!         
+            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_s, zua_ups, zva_ups )
+            !
+            zamsk = 1._wp
+            !== Salt content ==!
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+               &                                      psv_i, psv_i )
+            !== Ice heat content ==!
+            DO jk = 1, nlay_i
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+                  &                                      pe_i(:,:,jk,:), pe_i(:,:,jk,:) )
+            END DO
+            !== Snw heat content ==!
+            DO jk = 1, nlay_s
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+                  &                                      pe_s(:,:,jk,:), pe_s(:,:,jk,:) )
+            END DO
+            !
+            !                          ! ------------------------------------------ !
+         ELSEIF( np_advS == 2 ) THEN   ! -- advection form: -div( uA * uHS / u ) -- !
+            !                          ! ------------------------------------------ !
+            zamsk = 0._wp
+            !== Ice volume ==!
+            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_i, zua_ups, zva_ups )
+            !== Snw volume ==!         
+            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_s, zua_ups, zva_ups )
+            !== Salt content ==!
+            zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, psv_i, zua_ups, zva_ups )
+            !== Ice heat content ==!
+            DO jk = 1, nlay_i
+               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups )
+            END DO
+            !== Snw heat content ==!
+            DO jk = 1, nlay_s
+               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups )
+            END DO
+            !
+            !                          ! ----------------------------------------- !
+         ELSEIF( np_advS == 3 ) THEN   ! -- advection form: -div( uV * uS / u ) -- !
+            !                          ! ----------------------------------------- !
+            zamsk = 0._wp
+            !
+            ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl),  &
+               &      zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) )
+            !
+            ! inverse of Vi
+            WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:)
+            ELSEWHERE                        ;   z1_vi(:,:,:) = 0.
+            END WHERE
+            ! inverse of Vs
+            WHERE( pv_s(:,:,:) >= epsi20 )   ;   z1_vs(:,:,:) = 1._wp / pv_s(:,:,:)
+            ELSEWHERE                        ;   z1_vs(:,:,:) = 0.
+            END WHERE
+            !
+            ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays)
+            !
+            !== Ice volume ==!
+            zuv_ups = zua_ups
+            zvv_ups = zva_ups
+            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
+            !== Salt content ==!
+            zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
+               &                                      zhvar, psv_i, zuv_ups, zvv_ups )
+            !== Ice heat content ==!
+            DO jk = 1, nlay_i
+               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups )
+            END DO
+            !== Snow volume ==!         
+            zuv_ups = zua_ups
+            zvv_ups = zva_ups
+            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
+            !== Snw heat content ==!
+            DO jk = 1, nlay_s
+               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
+                  &                                      zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups )
+            END DO
+            !
+            DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs )
+            !
+         ENDIF
+         !
+         !== Ice age ==!
+         zamsk = 1._wp
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+            &                                      poa_i, poa_i )
+         !
+         !== melt ponds ==!
          IF ( ln_pnd_H12 ) THEN
-            ! set u*a=u for advection of Ap only 
-            DO jl = 1, jpl
-               zua_ho(:,:,jl) = zudy(:,:)
-               zva_ho(:,:,jl) = zvdx(:,:)
-            END DO
-            
+            ! concentration
             zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, &
+               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
+            ! volume
             zamsk = 0._wp
-            !
             zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! mp volume
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_ip, zua_ups, zva_ups )
          ENDIF
          !
+         !== Open water area ==!
          zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                   ! Open water area
+               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 
                   &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
             END DO
          END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. )
-         !
+         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. )
+         !
+         !
+         ! --- Ensure non-negative fields and in-bound thicknesses --- !
+         ! 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( 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( 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
       !
@@ -214,5 +361,6 @@
 
    
-   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
+   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox,  &
+      &                                            pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -221,40 +369,51 @@
       !!                 tracers and add it to the general trend of tracer equations
       !!
-      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracer H
+      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc
       !!                 - calculate tracer H at u and v points (Ultimate)
-      !!                 - calculate the high order fluxes using alterning directions (Macho?)
+      !!                 - calculate the high order fluxes using alterning directions (Macho)
       !!                 - apply a limiter on the fluxes (nonosc_ice)
-      !!                 - convert this tracer flux to a tracer content flux (uH -> uV)
-      !!                 - calculate the high order solution for tracer content V
+      !!                 - convert this tracer flux to a "volume" flux (uH -> uV)
+      !!                 - apply a limiter a second time on the volumes fluxes (nonosc_ice)
+      !!                 - calculate the high order solution for V
       !!
-      !! ** Action : solve 2 equations => a) da/dt = -div(ua)
-      !!                                  b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
-      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H.
-      !!                        - then we convert this flux to a "volume" flux this way => uH*ua/u
-      !!                             where ua is the flux from eq. a)
-      !!                        - at last we estimate dV/dt = -div(uH*ua/u)
+      !! ** Action : solve 3 equations => a) dA/dt  = -div(uA)
+      !!                                  b) dV/dt  = -div(uV)  using dH/dt = -u.grad(H)
+      !!                                  c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S)
       !!
-      !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc.
-      !!           - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now
-      !!             is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of
-      !!             the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
+      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H.
+      !!                        - then we convert this flux to a "volume" flux this way => uH * uA / u
+      !!                             where uA is the flux from eq. a)
+      !!                             this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur)
+      !!                        - at last we estimate dV/dt = -div(uH * uA / u)
+      !!
+      !!             in eq. c), one can solve the equation for  S (ln_advS=T), then dVS/dt = -div(uV * uS  / u)
+      !!                                                or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u) 
+      !!
+      !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc.
+      !!           - At the ice edge, Ultimate scheme can lead to:
+      !!                              1) negative interpolated tracers at u-v points
+      !!                              2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward
+      !!                              Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of
+      !!                                               the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
+      !!                              Solution for 2): we set it to 0 in this case
       !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
       !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
-      !!             work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D.
+      !!             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).
-      !! To-do: expand the prelimiter from zalesak to make it work in 2D
-      !!----------------------------------------------------------------------
-      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)
-      INTEGER                         , INTENT(in   )           ::   jt             ! number of sub-iteration
-      INTEGER                         , INTENT(in   )           ::   kt             ! number of iteration
-      REAL(wp)                        , INTENT(in   )           ::   pdt            ! tracer time-step
-      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
-      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
-      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
-      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt             ! tracer field
-      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc            ! tracer content field
-      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
+      !!             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.
+      !!----------------------------------------------------------------------
+      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)
+      INTEGER                         , INTENT(in   )           ::   jt               ! number of sub-iteration
+      INTEGER                         , INTENT(in   )           ::   kt               ! number of iteration
+      REAL(wp)                        , INTENT(in   )           ::   pdt              ! tracer time-step
+      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv        ! 2 ice velocity components => u*e2
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc       ! 2 ice velocity components => u*e2 or u*a*e2u
+      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox     ! upstream velocity
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt               ! tracer field
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc              ! tracer content field
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL ::   pua_ups, pva_ups ! upstream u*a fluxes
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho   ! high order u*a fluxes
       !
       INTEGER  ::   ji, jj, jl       ! dummy loop indices  
@@ -289,7 +448,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
-                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
-                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
+                  IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
+                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj)
+                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
                   ELSE
                      zfu_ho (ji,jj,jl) = 0._wp
@@ -297,7 +456,7 @@
                   ENDIF
                   !
-                  IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
-                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
-                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
+                  IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
+                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
+                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
                   ELSE
                      zfv_ho (ji,jj,jl) = 0._wp  
@@ -307,15 +466,36 @@
             END DO
          END DO
+
+         ! the new "volume" fluxes must also be "flux corrected"
+         ! thus we calculate the upstream solution and apply a limiter again
+         DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
+                  !
+                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. )
+         !
+         IF    ( np_limiter == 1 ) THEN
+            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
+         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
+            CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho )
+            CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho )
+         ENDIF
+         !
       ENDIF
-      !                                   --ho
-      ! in case of advection of A: output u*a
-      ! -------------------------------------
+      !                                   --ho    --ups
+      ! in case of advection of A: output u*a and u*a
+      ! -----------------------------------------------
       IF( PRESENT( pua_ho ) ) THEN
          DO jl = 1, jpl
             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)
-               END DO
+                  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)
+              END DO
             END DO
          END DO
@@ -485,7 +665,7 @@
          END DO
          !
-         IF    ( kn_limiter == 1 ) THEN
+         IF    ( np_limiter == 1 ) THEN
             CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
-         ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
+         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
             CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
             CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
@@ -503,5 +683,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
@@ -524,5 +704,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
 
          ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
@@ -535,5 +715,5 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
@@ -556,8 +736,8 @@
                END DO
             END DO
-            IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
 
          ENDIF
-         IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
          
       ENDIF
@@ -595,7 +775,7 @@
          !
          !                                                        !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )
+         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
          !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+         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
@@ -605,5 +785,5 @@
                      &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) &
                      &                                                                                        * pamsk           &
-                     &                             ) * pdt ) * tmask(ji,jj,1) 
+                     &                             ) * pdt ) * tmask(ji,jj,1)
                END DO
             END DO
@@ -613,10 +793,10 @@
          !                                                        !--  ultimate interpolation of pt at v-point  --!
          IF( ll_hoxy ) THEN
-            CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
+            CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
          ELSE
-            CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )
+            CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
          ENDIF
          !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
          !         
          !
@@ -624,7 +804,7 @@
          !
          !                                                        !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )
+         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
          !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+         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
@@ -642,19 +822,19 @@
          !                                                        !--  ultimate interpolation of pt at u-point  --!
          IF( ll_hoxy ) THEN
-            CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
+            CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
          ELSE
-            CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )
+            CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
          ENDIF
          !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
          !
       ENDIF
 
-      IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+      IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
       !
    END SUBROUTINE macho
 
 
-   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )
+   SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -666,4 +846,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
+      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)
       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
@@ -792,5 +973,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN
+                  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)   &
                         &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -812,5 +993,5 @@
    
  
-   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )
+   SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_y  ***
@@ -822,4 +1003,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
+      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)
       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
@@ -945,5 +1127,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN
+                  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) )  &
                         &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1009,5 +1191,5 @@
       !                        |      |      |        |    *
       !            t_ups :       i-1     i       i+1       i+2   
-      IF( ll_prelimiter_zalesak ) THEN
+      IF( ll_prelim ) THEN
          
          DO jl = 1, jpl
@@ -1088,10 +1270,11 @@
                !
                !                                  ! up & down beta terms
-               IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
-               ELSE                    ; zbetup(ji,jj,jl) = 0._wp ! zbig
+               ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
+               IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
+               ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig
                ENDIF
                !
-               IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
-               ELSE                    ; zbetdo(ji,jj,jl) = 0._wp ! zbig
+               IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
+               ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig
                ENDIF
                !
@@ -1135,18 +1318,4 @@
          END DO
 
-         ! clem test
-!!         DO jj = 2, jpjm1
-!!            DO ji = 2, fs_jpim1   ! vector opt.
-!!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
-!!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
-!!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
-!!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
-!!                  &         ) * tmask(ji,jj,1)
-!!               IF( zzt < -epsi20 ) THEN
-!!                  WRITE(numout,*) 'T<0 nonosc_ice',zzt
-!!               ENDIF
-!!            END DO
-!!         END DO
-
       END DO
       !
@@ -1189,5 +1358,5 @@
                Rjp = zslpx(ji+1,jj,jl)
 
-               IF( kn_limiter == 3 ) THEN
+               IF( np_limiter == 3 ) THEN
 
                   IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
@@ -1205,5 +1374,5 @@
                   pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
 
-               ELSEIF( kn_limiter == 2 ) THEN
+               ELSEIF( np_limiter == 2 ) THEN
                   IF( Rj /= 0. ) THEN
                      IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
@@ -1284,5 +1453,5 @@
                Rjp = zslpy(ji,jj+1,jl)
 
-               IF( kn_limiter == 3 ) THEN
+               IF( np_limiter == 3 ) THEN
 
                   IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
@@ -1300,5 +1469,5 @@
                   pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
 
-               ELSEIF( kn_limiter == 2 ) THEN
+               ELSEIF( np_limiter == 2 ) THEN
 
                   IF( Rj /= 0. ) THEN
@@ -1344,4 +1513,94 @@
    END SUBROUTINE limiter_y
 
+
+   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  ***
+      !!
+      !! ** Purpose : Thickness correction in case advection scheme creates
+      !!              abnormally tick ice or snow
+      !!
+      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by adapting ice concentration
+      !!              2- check whether snow thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by sending the excess in the ocean
+      !!              3- check whether snow load deplets the snow-ice interface below sea level$
+      !!                 and reduce it by sending the excess in the ocean
+      !!              4- correct pond concentration to avoid a_ip > a_i
+      !!
+      !! ** 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
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i
+      !
+      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
+      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
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+                  !
+                  !                               ! -- check h_ip -- !
+                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
+                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
+                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
+                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
+                     ENDIF
+                  ENDIF
+                  !
+                  !                               ! -- check h_i -- !
+                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
+                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
+                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
+                  ENDIF
+                  !
+                  !                               ! -- check h_s -- !
+                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
+                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
+                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+                     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 * 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
+                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
+                  ENDIF           
+                  !
+                  !                               ! -- check snow load -- !
+                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
+                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
+                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
+                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+                  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 * 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
+                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
+                  ENDIF
+                  
+               ENDIF
+            END DO
+         END DO
+      END DO 
+      !                                           !-- correct pond concentration to avoid a_ip > a_i
+      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
+      !
+      !
+   END SUBROUTINE Hbig
+   
 #else
    !!----------------------------------------------------------------------
