Index: /NEMO/releases/release-4.0/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml
===================================================================
--- /NEMO/releases/release-4.0/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 10909)
+++ /NEMO/releases/release-4.0/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 10910)
@@ -79,9 +79,4 @@
        <field field_ref="vfxsnw"           name="vfxsnw" />
        
-       <!-- diag error for negative ice volume after advection -->
-       <field field_ref="iceneg_pres"      name="sineg_pres" />
-       <field field_ref="iceneg_volu"      name="sineg_volu" />
-       <field field_ref="iceneg_hfx"       name="sineg_hfx"  />
-
        <!-- categories -->
        <field field_ref="icemask_cat"      name="simskcat"/>
Index: /NEMO/releases/release-4.0/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/releases/release-4.0/cfgs/SHARED/field_def_nemo-ice.xml	(revision 10909)
+++ /NEMO/releases/release-4.0/cfgs/SHARED/field_def_nemo-ice.xml	(revision 10910)
@@ -160,9 +160,4 @@
           <field id="hfxdhc"       long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" />
 
-	  <!-- diagnostics of the negative values resulting from the advection scheme -->
-	  <field id="iceneg_pres"  long_name="Fraction of time steps with negative sea ice volume"                    unit=""     />
-	  <field id="iceneg_volu"  long_name="Negative sea ice volume per area arising from advection"                unit="m"    />
-          <field id="iceneg_hfx"   long_name="Negative sea ice heat content (eq. heat flux) arising from advection"   unit="W/m2" />
-
 	  <!-- sbcssm variables -->
           <field id="sst_m"    unit="degC" />
Index: /NEMO/releases/release-4.0/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/releases/release-4.0/cfgs/SHARED/namelist_ice_ref	(revision 10909)
+++ /NEMO/releases/release-4.0/cfgs/SHARED/namelist_ice_ref	(revision 10910)
@@ -154,5 +154,5 @@
 &namthd_do      !   Ice growth in open water
 !------------------------------------------------------------------------------
-   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice
+   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_himin
    ln_frazil        = .false.         !  Frazil ice parameterization (ice collection as a function of wind)
       rn_maxfraz    =   1.0           !     maximum fraction of frazil ice collecting at the ice base
Index: /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml
===================================================================
--- /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 10909)
+++ /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 10910)
@@ -79,9 +79,4 @@
        <field field_ref="vfxsnw"           name="vfxsnw" />
 
-       <!-- diag error for negative ice volume after advection -->
-       <field field_ref="iceneg_pres"      name="sineg_pres" />
-       <field field_ref="iceneg_volu"      name="sineg_volu" />
-       <field field_ref="iceneg_hfx"       name="sineg_hfx"  />
-       
        <!-- categories -->
        <field field_ref="icemask_cat"      name="simskcat"/>
Index: /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_cfg
===================================================================
--- /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 10909)
+++ /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 10910)
@@ -80,4 +80,6 @@
                      ! Sea-ice :
    nn_ice      = 2         !  SI3
+   ln_ice_embd = .false.   !  =T embedded sea-ice (pressure + mass and salt exchanges)
+      !                    !  =F levitating ice (no pressure, mass and salt exchanges)
                      ! Misc. options of sbc :
    ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr )
@@ -90,5 +92,5 @@
       !
       ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012)
-      ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015)
+      ln_Cd_L15   = .true.    !  air-ice drags = F(ice concentration) (Lupkes et al. 2015)
       !
    cn_dir = './'  !  root directory for the bulk data location
@@ -96,16 +98,4 @@
    !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !          weights filename            ! rotation ! land/sea mask !
    !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   !
-!! ERAI
-!   sn_wndi     = 'u10_era_spitz'               ,         6         , 'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , ''
-!   sn_wndj     = 'v10_era_spitz'               ,         6         , 'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , ''
-!   sn_qsr      = 'ssrd_era_spitz'              ,         6         , 'ssrd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_qlw      = 'strd_era_spitz'              ,         6         , 'strd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_tair     = 't2_era_spitz'                ,         6         , 't2'      ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_humi     = 'humi_era_spitz'              ,         6         , 'humi'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_prec     = 'precip_era_spitz'            ,         6         , 'precip'  ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_snow     = 'snow_era_spitz'              ,         6         , 'snow'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , ''
-!   sn_tdif     = 'taudif'                      ,         6         , 'taudif'  ,   .true.    , .false. , 'yearly'  , ''             , '' , ''
-!      rn_zqt      = 10.       !  Air temperature & humidity reference height (m)
-!! MAR
    sn_wndi     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , ''
    sn_wndj     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , ''
@@ -129,12 +119,4 @@
    !                       ! type of penetration                        (default: NO selection)
    ln_qsr_rgb  = .true.       !  RGB light penetration (Red-Green-Blue)
-/
-!-----------------------------------------------------------------------
-&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T)
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namsbc_rnf    !   runoffs                                              (ln_rnf =T)
-!-----------------------------------------------------------------------
 /
 !!======================================================================
@@ -236,19 +218,13 @@
 !-----------------------------------------------------------------------
    ln_loglayer = .true.   !  logarithmic drag: Cd = vkarmn/log(z/z0) |U|
+   ln_drgimp   = .true.   !  implicit top/bottom friction flag
 /
 !-----------------------------------------------------------------------
 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F)
 !-----------------------------------------------------------------------
-   rn_Cd0      =  2.5e-3   !!1.e-3    !  drag coefficient [-]
-   rn_Cdmax    =  0.01     !!0.1      !  drag value maximum [-] (logarithmic drag)
-   rn_ke0      =  0.       !!2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases)
+   rn_Cd0      =  2.5e-3   !  drag coefficient [-]
+   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag)
+   rn_ke0      =  0.       !  background kinetic energy  [m2/s2] (non-linear cases)
    rn_z0       =  3.e-3    !  roughness [m] (ln_loglayer=T)
-   ln_boost    = .false.   !  =T regional boost of Cd0 ; =F constant
-      rn_boost =  50.         !  local boost factor  [-]
-/
-!-----------------------------------------------------------------------
-&nambbc        !   bottom temperature boundary condition                (default: OFF)
-!-----------------------------------------------------------------------
-   ln_trabbc   = .false.    !  Apply a geothermal heating at the ocean bottom
 /
 !-----------------------------------------------------------------------
@@ -258,6 +234,4 @@
       nn_bbl_ldf  =  1        !  diffusive bbl (=1)   or not (=0)
       nn_bbl_adv  =  0        !  advective bbl (=1/2) or not (=0)
-      rn_ahtbbl   =  1000.    !  lateral mixing coefficient in the bbl  [m2/s]
-      rn_gambbl   =  10.      !  advective bbl coefficient                 [s]
 /
 !!======================================================================
@@ -293,5 +267,5 @@
    !                       !  Coefficients:
    nn_aht_ijk_t    = 31        !  space/time variation of eddy coefficient:
-      !                             !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing)
+   !                                !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing)
 /
 !!======================================================================
@@ -320,5 +294,5 @@
 &namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynvor_enT = .true. !  energy conserving scheme (T-point)
+   ln_dynvor_eeT = .true.  !  energy conserving scheme (een using e3t)
 /
 !-----------------------------------------------------------------------
@@ -352,8 +326,14 @@
 &namzdf        !   vertical physics manager                             (default: NO selection)
 !-----------------------------------------------------------------------
+   !                       ! adaptive-implicit vertical advection
+   ln_zad_Aimp = .true.      !  Courant number dependent scheme (Shchepetkin 2015)
    !                       ! type of vertical closure (required)
    ln_zdftke   = .true.       !  Turbulent Kinetic Energy closure       (T =>   fill namzdf_tke)
    ln_zdfgls   = .false.      !  Generic Length Scale closure           (T =>   fill namzdf_gls)
+   !                       ! convection
+   ln_zdfevd   = .true.       !  enhanced vertical diffusion
+   !
    ln_zdfddm   = .true.    ! double diffusive mixing
+   !
    !                       !  Coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
@@ -384,8 +364,4 @@
 /
 !-----------------------------------------------------------------------
-&nam_diaharm   !   Harmonic analysis of tidal constituents              ('key_diaharm')
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4")
 !-----------------------------------------------------------------------
Index: /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_ice_cfg
===================================================================
--- /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_ice_cfg	(revision 10909)
+++ /NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_ice_cfg	(revision 10910)
@@ -35,4 +35,5 @@
 &namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
+   ln_landfast_L16  = .true.          !  landfast: parameterization from Lemieux 2016
 /
 !------------------------------------------------------------------------------
@@ -43,4 +44,6 @@
 &namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
+   ln_rhg_EVP       = .true.          !  EVP rheology
+      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017)
 /
 !------------------------------------------------------------------------------
@@ -67,5 +70,5 @@
 &namthd_do      !   Ice growth in open water
 !------------------------------------------------------------------------------
-   rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_hnewice
+   rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_himin
    ln_frazil        = .true.          !  Frazil ice parameterization (ice collection as a function of wind)
 /
@@ -77,7 +80,6 @@
 &namthd_pnd     !   Melt ponds
 !------------------------------------------------------------------------------
-   ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
-   ln_pnd_CST       = .false.         !  activate constant melt ponds
-   ln_pnd_alb       = .false.         !  melt ponds affect albedo or not
+   ln_pnd_H12       = .true.          !  activate evolutive melt ponds (from Holland et al 2012)
+   ln_pnd_alb       = .true.          !  melt ponds affect albedo or not
 /
 
Index: /NEMO/releases/release-4.0/src/ICE/icedyn.F90
===================================================================
--- /NEMO/releases/release-4.0/src/ICE/icedyn.F90	(revision 10909)
+++ /NEMO/releases/release-4.0/src/ICE/icedyn.F90	(revision 10910)
@@ -74,8 +74,7 @@
       INTEGER, INTENT(in) ::   kt     ! ice time step
       !!
-      INTEGER  ::   ji, jj, jl        ! dummy loop indices
+      INTEGER  ::   ji, jj        ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv
-      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i
       !!--------------------------------------------------------------------
       !
@@ -89,48 +88,37 @@
       ENDIF
       !                      
-      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
-         tau_icebfr(:,:) = 0._wp
-         DO jl = 1, jpl
-            WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
-         END DO
-      ENDIF
-      !
-      !                                !-- 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, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), &
-                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), &
-                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), &
-                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) )
-               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), &
-                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), &
-                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), &
-                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) )
-               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), &
-                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), &
-                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), &
-                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
-      !
-      !
-      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
+      ! retrieve thickness from volume for landfast param. and UMx advection scheme
+      WHERE( a_i(:,:,:) >= epsi20 )
+         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
+         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
+      ELSEWHERE
+         h_i(:,:,:) = 0._wp
+         h_s(:,:,:) = 0._wp
+      END WHERE
+      !
+      WHERE( a_ip(:,:,:) >= epsi20 )
+         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
+      ELSEWHERE
+         h_ip(:,:,:) = 0._wp
+      END WHERE
+      !
+      !
+      SELECT CASE( nice_dyn )          !-- Set which dynamics is running
 
       CASE ( np_dynALL )           !==  all dynamical processes  ==!
-         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
-         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting 
-         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections
-
+         !
+         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology  
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting 
+         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
+         !
       CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
-         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
-         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting)
-
+         !
+         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology  
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
+         !
       CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
-         ALLOCATE( zdivu_i(jpi,jpj) )
+         !
          ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
          ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
@@ -145,20 +133,8 @@
          END DO
          ! ---
-         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
-
-         ! diagnostics: divergence at T points 
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
-         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
-
-         DEALLOCATE( zdivu_i )
-
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+         !
       CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
-         ALLOCATE( zdivu_i(jpi,jpj) )
+         !
          u_ice(:,:) = rn_uice * umask(:,:,1)
          v_ice(:,:) = rn_vice * vmask(:,:,1)
@@ -166,19 +142,30 @@
          !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
          ! ---
-         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
-
-         ! diagnostics: divergence at T points 
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
+
+      END SELECT
+      !
+      !
+      ! diagnostics: divergence at T points 
+      IF( iom_use('icediv') ) THEN
+         !
+         SELECT CASE( nice_dyn )
+
+         CASE ( np_dynADV1D , np_dynADV2D )
+
+            ALLOCATE( zdivu_i(jpi,jpj) )
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+               END DO
             END DO
-         END DO
-         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
-         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
-
-         DEALLOCATE( zdivu_i )
-
-      END SELECT
+            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
+            CALL iom_put( "icediv" , zdivu_i(:,:) )
+            DEALLOCATE( zdivu_i )
+
+         END SELECT
+         !
+      ENDIF
       !
       ! controls
@@ -186,92 +173,4 @@
       !
    END SUBROUTINE ice_dyn
-
-
-   SUBROUTINE Hbig( phi_max, phs_max, phip_max )
-      !!-------------------------------------------------------------------
-      !!                  ***  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 fraction to avoid a_ip > a_i
-      !!
-      !! ** input   : Max thickness of the surrounding 9-points
-      !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
-      !
-      INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra
-      !!-------------------------------------------------------------------
-      ! controls
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-      CALL ice_var_zapsmall                       !-- zap small areas
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( v_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. v_ip(ji,jj,jl) > 0._wp ) THEN
-                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )
-                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN
-                        a_ip(ji,jj,jl) = v_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 = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
-                     a_i(ji,jj,jl) = v_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 = v_s(ji,jj,jl) / a_i(ji,jj,jl)
-                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
-                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
-                     !
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
-                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
-                     !
-                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
-                     v_s(ji,jj,jl)          = a_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, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
-                  IF( zvs_excess > 0._wp ) THEN
-                     zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_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( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
-                     !
-                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
-                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess
-                  ENDIF
-                  
-               ENDIF
-            END DO
-         END DO
-      END DO 
-      !                                           !-- correct pond fraction to avoid a_ip > a_i
-      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
-      !
-      ! controls
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      !
-   END SUBROUTINE Hbig
 
 
@@ -337,16 +236,16 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn:'
-         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL
-         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV
-         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D
-         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D
-         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')'
-         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat
-         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16
-         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home
-         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra
-         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr
-         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax
-         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile
+         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL
+         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV
+         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D
+         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D
+         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')'
+         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
+         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
+         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home
+         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
+         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
+         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
+         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
          WRITE(numout,*)
       ENDIF
Index: /NEMO/releases/release-4.0/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/releases/release-4.0/src/ICE/icedyn_adv.F90	(revision 10909)
+++ /NEMO/releases/release-4.0/src/ICE/icedyn_adv.F90	(revision 10910)
@@ -64,10 +64,9 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! number of iteration
-      !
-      INTEGER ::   jl   ! dummy loop indice
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmask  ! fraction of time step with v_i < 0
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('icedyn_adv')
+      ! controls
+      IF( ln_timing    )   CALL timing_start('icedyn_adv')                                                             ! timing
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -76,44 +75,20 @@
          WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
-      
-      ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-                     
-      !----------
-      ! Advection
-      !----------
+      !
+      !---------------!
+      !== Advection ==!
+      !---------------!
       SELECT CASE( nice_adv )
       !                                !-----------------------!
       CASE( np_advUMx )                ! ULTIMATE-MACHO scheme !
          !                             !-----------------------!
-         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
-      !                                !-----------------------!
+         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, &
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+         !                             !-----------------------!
       CASE( np_advPRA )                ! PRATHER scheme        !
          !                             !-----------------------!
-         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, &
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
       END SELECT
-
-      !----------------------------
-      ! Debug the advection schemes
-      !----------------------------
-      ! clem: At least one advection scheme above is not strictly positive => UMx
-      !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes)
-      !       In UMx    , advected fields are not perfectly bounded and negative values can appear.
-      !                   These values are usually very small but in some occasions they can also be non-negligible
-      !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix)
-      !
-      ! record the negative values resulting from UMx
-      zmask(:,:) = 0._wp ! keep the init to 0 here
-      DO jl = 1, jpl
-         WHERE( v_i(:,:,jl) < 0._wp )   zmask(:,:) = 1._wp
-      END DO
-      IF( iom_use('iceneg_pres') )   CALL iom_put("iceneg_pres", zmask                                      )  ! fraction of time step with v_i < 0
-      IF( iom_use('iceneg_volu') )   CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 )                )  ! negative ice volume (only)
-      IF( iom_use('iceneg_hfx' ) )   CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. )  &  ! negative ice heat content (only)
-         &                                                                  , dim=4 ), dim=3 ) )* r1_rdtice )  ! -- eq. heat flux [W/m2]
-      !
-      ! ==> conservation is ensured by calling this routine below,
-      !     however the global ice volume is then changed by advection (but errors are small) 
-      CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
 
       !------------
Index: /NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90	(revision 10909)
+++ /NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90	(revision 10910)
@@ -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
@@ -39,7 +39,11 @@
    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 
+   ! 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.
+
+   ! 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.
    
@@ -51,7 +55,9 @@
    
    ! 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
 
+   ! 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
@@ -64,5 +70,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 )
       !!----------------------------------------------------------------------
@@ -79,4 +85,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
@@ -94,17 +103,41 @@
       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(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) ::   zua_ho, zva_ho, zua_ups, zva_ups
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups 
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip, z1_vi, z1_vs 
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
       !!----------------------------------------------------------------------
       !
       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(:,:) ) )
@@ -116,5 +149,4 @@
       ELSE                         ;   icycle = 1
       ENDIF
-      
       zdt = rdt_ice / REAL(icycle)
 
@@ -154,34 +186,95 @@
          END WHERE
          !
-         ! set u*a=u for advection of A only 
+         ! set u*A=u for advection of A only 
          DO jl = 1, jpl
             zua_ho(:,:,jl) = zudy(:,:)
             zva_ho(:,:,jl) = zvdx(:,:)
          END DO
-         
+
+         !== 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
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
          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
-         !
-         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
-         !
-         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN                                                              !-- Ice Age
+
+         IF( .NOT. ll_advS ) THEN   !-- advection form: uA * uHS / u
+            !== 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
+            !
+         ELSE                       !-- advection form: uV * uS / u
+            ! 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
+            !
+         ENDIF
+         !
+         !== Ice age ==!
+         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
             ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with
             !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that
@@ -189,5 +282,5 @@
             !!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 )
-            ! set u*a=u for advection of ice age
+            ! set u*A=u for advection of ice age
             DO jl = 1, jpl
                zua_ho(:,:,jl) = zudy(:,:)
@@ -195,27 +288,32 @@
             END DO
             zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i )
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho, zva_ho, zcu_box, zcv_box, &
+               &                                      poa_i, poa_i )
             zamsk = 0._wp
          ENDIF
          !
-         IF ( ln_pnd_H12 ) THEN                                                                                               !-- melt ponds
-            ! set u*a=u for advection of Ap only 
+         !== 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
-            !
+            ! fraction
             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 ) ! fraction
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
             zamsk = 0._wp
-            !
+            ! volume
             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 )                 ! 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
@@ -223,4 +321,15 @@
          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( 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 )
+
       END DO
       !
@@ -228,5 +337,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  ***
@@ -235,40 +345,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  
@@ -303,7 +424,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
@@ -311,7 +432,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  
@@ -321,15 +442,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    ( kn_limiter == 1 ) THEN
+            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
+         ELSEIF( kn_limiter == 2 .OR. kn_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
@@ -609,5 +751,5 @@
          !
          !                                                        !--  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 )
@@ -619,5 +761,8 @@
                      &                              + 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)
+                  !!clem test
+                  !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) )
+                  !!clem test
                END DO
             END DO
@@ -627,7 +772,7 @@
          !                                                        !--  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 --!
@@ -638,5 +783,5 @@
          !
          !                                                        !--  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 )
@@ -649,4 +794,7 @@
                      &                                                                                        * 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
@@ -656,7 +804,7 @@
          !                                                        !--  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 --!
@@ -670,5 +818,5 @@
 
 
-   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  ***
@@ -680,4 +828,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
@@ -688,5 +837,5 @@
       !
       INTEGER  ::   ji, jj, jl             ! dummy loop indices
-      REAL(wp) ::   zcu, zdx2, zdx4        !   -      -
+      REAL(wp) ::   zcu, zdx2, zdx4, zvi_cen2        !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4
       !!----------------------------------------------------------------------
@@ -799,4 +948,19 @@
       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
@@ -806,5 +970,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN
+                  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
                      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) ) )
@@ -826,5 +991,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  ***
@@ -836,4 +1001,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
@@ -844,5 +1010,5 @@
       !
       INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
+      REAL(wp) ::   zcv, zdy2, zdy4, zvi_cen2    !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4
       !!----------------------------------------------------------------------
@@ -952,4 +1118,19 @@
       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
@@ -959,5 +1140,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN
+                  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
                      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) ) )
@@ -1102,10 +1284,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
                !
@@ -1148,18 +1331,4 @@
             END DO
          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
@@ -1358,4 +1527,92 @@
    END SUBROUTINE limiter_y
 
+
+   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 )
+      !!-------------------------------------------------------------------
+      !!                  ***  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 fraction to avoid a_ip > a_i
+      !!
+      !! ** input   : Max thickness of the surrounding 9-points
+      !!-------------------------------------------------------------------
+      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) ::   zhip, zhi, zhs, zvs_excess, zfra
+      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
+      !!-------------------------------------------------------------------
+      !
+      !
+      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 * 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
+                     !
+                     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 * 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
+                     !
+                     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 fraction to avoid a_ip > a_i
+      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
+      !
+      !
+   END SUBROUTINE Hbig
+   
 #else
    !!----------------------------------------------------------------------
Index: /NEMO/releases/release-4.0/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/releases/release-4.0/src/ICE/icedyn_rhg.F90	(revision 10909)
+++ /NEMO/releases/release-4.0/src/ICE/icedyn_rhg.F90	(revision 10910)
@@ -58,4 +58,6 @@
       !!--------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ice time step
+      !
+      INTEGER  ::   jl   ! dummy loop indices
       !!--------------------------------------------------------------------
       ! controls
@@ -68,8 +70,15 @@
          WRITE(numout,*)'~~~~~~~~~~~'
       ENDIF
-
-      ! --------
-      ! Rheology
-      ! --------   
+      !
+      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
+         tau_icebfr(:,:) = 0._wp
+         DO jl = 1, jpl
+            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
+         END DO
+      ENDIF
+      !
+      !--------------!
+      !== Rheology ==!
+      !--------------!   
       SELECT CASE( nice_rhg )
       !                                !------------------------!
Index: /NEMO/releases/release-4.0/src/ICE/icewri.F90
===================================================================
--- /NEMO/releases/release-4.0/src/ICE/icewri.F90	(revision 10909)
+++ /NEMO/releases/release-4.0/src/ICE/icewri.F90	(revision 10910)
@@ -145,5 +145,5 @@
 
         ! record presence of fast ice
-         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
+         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
          ELSEWHERE                                                ; zfast(:,:) = 0._wp
          END WHERE
