Index: NEMO/trunk/src/ICE/ice.F90
===================================================================
--- NEMO/trunk/src/ICE/ice.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/ice.F90	(revision 10413)
@@ -129,10 +129,15 @@
    !                                     !!** ice-dynamics namelist (namdyn) **
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
-   LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F) 
-   REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice
-   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (landfast ice) 
-   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice) 
-   !
-   !                                     !!** ice-rheology namelist (namrhg) **
+   LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016 
+   LOGICAL , PUBLIC ::   ln_landfast_home !: landfast ice parameterizationfrom home made 
+   REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice
+   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home) 
+   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction
+   REAL(wp), PUBLIC ::   rn_tensile       !:    isotropic tensile strength
+   !
+   !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) **
+   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength (also used for landfast param)
+   !
+   !                                     !!** ice-rheology namelist (namdyn_rhg) **
    LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
@@ -140,4 +145,8 @@
    INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
+   !
+   !                                     !!** ice-advection namelist (namdyn_adv) **
+   LOGICAL , PUBLIC ::   ln_adv_Pra       !: Prather        advection scheme
+   LOGICAL , PUBLIC ::   ln_adv_UMx       !: Ultimate-Macho advection scheme
    !
    !                                     !!** ice-surface forcing namelist (namforcing) **
Index: NEMO/trunk/src/ICE/icedyn.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icedyn.F90	(revision 10413)
@@ -21,4 +21,5 @@
    USE icecor         ! sea-ice: corrections
    USE icevar         ! sea-ice: operations
+   USE icectl         ! sea-ice: control prints
    !
    USE in_out_manager ! I/O manager
@@ -37,14 +38,16 @@
    INTEGER ::              nice_dyn   ! choice of the type of dynamics
    !                                        ! associated indices:
-   INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   INTEGER, PARAMETER ::   np_dynALL     = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
    INTEGER, PARAMETER ::   np_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection) 
-   INTEGER, PARAMETER ::   np_dynADV     = 3   ! only advection w prescribed vel.(rn_uvice + advection)
+   INTEGER, PARAMETER ::   np_dynADV1D   = 3   ! only advection 1D - test case from Schar & Smolarkiewicz 1996
+   INTEGER, PARAMETER ::   np_dynADV2D   = 4   ! only advection 2D w prescribed vel.(rn_uvice + advection)
    !
    ! ** namelist (namdyn) **
-   LOGICAL  ::   ln_dynFULL       ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
-   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections  (rheology + advection)
-   LOGICAL  ::   ln_dynADV        ! only advection w prescribed vel.(rn_uvice + advection)
-   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV)
-   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV)
+   LOGICAL  ::   ln_dynALL        ! full ice dynamics                      (rheology + advection + ridging/rafting + correction)
+   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections         (rheology + advection)
+   LOGICAL  ::   ln_dynADV1D      ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996)
+   LOGICAL  ::   ln_dynADV2D      ! only advection in 2D w prescribed vel. (rn_uvice + advection)
+   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV1D & np_dynADV2D)
+   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV1D & np_dynADV2D)
    
    !! * Substitutions
@@ -53,5 +56,5 @@
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -71,9 +74,13 @@
       INTEGER, INTENT(in) ::   kt     ! ice time step
       !!
-      INTEGER ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax
+      INTEGER  ::   ji, jj, jl        ! dummy loop indices
+      REAL(wp) ::   zcoefu, zcoefv
+      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
       !!--------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('icedyn')
+      ! controls
+      IF( ln_timing    )   CALL timing_start('icedyn')                                                             ! timing
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -82,71 +89,122 @@
          WRITE(numout,*)'~~~~~~~'
       ENDIF
-
       !                      
-      IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction
+      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
          tau_icebfr(:,:) = 0._wp
          DO jl = 1, jpl
-            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_gamma )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
-         END DO
-         IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr )  
-      ENDIF
-
-      zhmax(:,:,:) = h_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
+            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
+               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( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1. )
+      !
+      !
+      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 )   ! -- advection of ice + correction on ice thickness
+         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 )   ! -- advection of ice + correction on ice thickness
+         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
+         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s 
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
+               zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
+               u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
+               v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
+            END DO
+         END DO
+         ! ---
+         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice + correction on ice thickness
+
+         ! diagnostics: divergence at T points 
+         DO jj = 2, jpjm1
             DO ji = 2, jpim1
-!!gm use of MAXVAL here is very probably less efficient than expending the 9 values
-               zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( h_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. )
-      !
-      !
-      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
-
-      CASE ( np_dynFULL )          !==  all dynamical processes  ==!
-         CALL ice_dyn_rhg   ( kt )                            ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness
-         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 )                            ! -- advection of ice
-         CALL Hpiling                                         ! -- simple pile-up (replaces ridging/rafting)
-
-      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities)
+               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( zdivu_i, 'T', 1. )
+         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
+
+         DEALLOCATE( zdivu_i )
+
+      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)
-         !!CALL RANDOM_NUMBER(u_ice(:,:))
-         !!CALL RANDOM_NUMBER(v_ice(:,:))
-         CALL ice_dyn_adv   ( kt )                            ! -- advection of ice
+         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
+         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
+         ! ---
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness
+
+         ! 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( zdivu_i, 'T', 1. )
+         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
+
+         DEALLOCATE( zdivu_i )
 
       END SELECT
-      !
-      IF( ln_timing )   CALL timing_stop('icedyn')
+       !
+      ! controls
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_timing    )   CALL timing_stop ('icedyn')                                                             ! timing
       !
    END SUBROUTINE ice_dyn
 
 
-   SUBROUTINE Hbig( phmax )
+   SUBROUTINE Hbig( phi_max, phs_max )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE Hbig  ***
       !!
       !! ** Purpose : Thickness correction in case advection scheme creates
-      !!              abnormally tick ice
-      !!
-      !! ** Method  : 1- check whether ice thickness resulting from advection is
-      !!                 larger than the surrounding 9-points before advection
-      !!                 and reduce it if a) divergence or b) convergence & at_i>0.8
-      !!              2- bound ice thickness with hi_max (99m)
+      !!              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) ::   phmax   ! max ice thick from surrounding 9-pts
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max   ! max ice thick from surrounding 9-pts
       !
       INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zh, zdv
+      REAL(wp) ::   zhi, zhs, zvs_excess, zfra
       !!-------------------------------------------------------------------
       !
@@ -156,14 +214,43 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax
+               IF ( v_i(ji,jj,jl) > 0._wp ) THEN
                   !
-                  zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl)
-                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)  
-                  !
-                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
-                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
-                     a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
+                  !                               ! -- 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)
+!!clem                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)  
+!!clem                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
+!!clem                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
+                  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 = a_i(ji,jj,jl) * phs_max(ji,jj,jl) / MAX( v_s(ji,jj,jl), 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) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     !
+                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,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 thickness very high (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 = 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) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     !
+                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra
+                     v_s(ji,jj,jl)   = v_s(ji,jj,jl) - zvs_excess
+                  ENDIF
+                  
                ENDIF
             END DO
@@ -215,7 +302,7 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
-         &             rn_ishlat  ,                                            &
-         &             ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax
+      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
+         &             rn_ishlat ,                                                           &
+         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
       !!-------------------------------------------------------------------
       !
@@ -233,13 +320,16 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn:'
-         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL
-         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV
-         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV
-         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 (T or F)                                ln_landfast  = ', ln_landfast
-         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_gamma     = ', rn_gamma
-         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,*) '      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
@@ -247,9 +337,11 @@
       ioptio = 0 
       !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
-      IF( ln_dynFULL   ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF
+      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF
       !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
       IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
-      !      !--- advection only with prescribed ice velocities (from namelist)
-      IF( ln_dynADV    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF
+      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
+      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF
+      !      !--- advection 2D only with prescribed ice velocities (from namelist)
+      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF
       !
       IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
@@ -261,6 +353,10 @@
       ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
       ENDIF
-      !                                      !--- NO Landfast ice : set to zero once for all
-      IF( .NOT.ln_landfast )   tau_icebfr(:,:) = 0._wp
+      !                                      !--- Landfast ice
+      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp
+      !
+      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN
+         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' )
+      ENDIF
       !
       CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
Index: NEMO/trunk/src/ICE/icedyn_adv.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_adv.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icedyn_adv.F90	(revision 10413)
@@ -40,7 +40,5 @@
    !
    ! ** namelist (namdyn_adv) **
-   LOGICAL ::   ln_adv_Pra   ! Prather        advection scheme
-   LOGICAL ::   ln_adv_UMx   ! Ultimate-Macho advection scheme
-   INTEGER ::   nn_UMx       ! order of the UMx advection scheme   
+   INTEGER         ::   nn_UMx       ! order of the UMx advection scheme   
    !
    !! * Substitution
@@ -49,5 +47,5 @@
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -89,11 +87,9 @@
       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, 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
 
@@ -101,7 +97,7 @@
       ! Debug the advection schemes
       !----------------------------
-      ! clem: At least one advection scheme above is not strictly positive => UM from 3d to 5th order
-      !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems not)
-      !       In UM3-5  , advected fields are not bounded and negative values can appear.
+      ! 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)
@@ -118,5 +114,5 @@
       !
       ! ==> conservation is ensured by calling this routine below,
-      !     however the global ice volume is then changed by advection (but errors are very small) 
+      !     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/trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 10413)
@@ -14,5 +14,5 @@
    !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
    !!   macho             : ???
-   !!   nonosc_2d         : compute monotonic tracer fluxes by a non-oscillatory algorithm 
+   !!   nonosc            : compute monotonic tracer fluxes by a non-oscillatory algorithm 
    !!----------------------------------------------------------------------
    USE phycst         ! physical constant
@@ -20,4 +20,5 @@
    USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition
    USE ice            ! sea-ice variables
+   USE icevar         ! sea-ice: operations
    !
    USE in_out_manager ! I/O manager
@@ -34,4 +35,52 @@
    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
 
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_ai, amaxu, amaxv
+   
+   LOGICAL ll_dens
+
+   ! advect H all the way (and get V=H*A at the end)
+   LOGICAL :: ll_thickness = .FALSE.
+   
+   ! look for 9 points around in nonosc limiter  
+   LOGICAL :: ll_9points = .FALSE.  ! false=better h?
+
+   ! use HgradU in nonosc limiter  
+   LOGICAL :: ll_HgradU = .TRUE.   ! no effect?
+
+   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
+   LOGICAL :: ll_neg = .TRUE.      ! keep TRUE
+   
+   ! limit the fluxes
+   LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !!
+   LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D
+   LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D
+   LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D
+
+   ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc)
+   LOGICAL :: ll_clem   = .TRUE.   ! simpler than rachid and works
+
+   ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*)
+   LOGICAL :: ll_gurvan = .FALSE.  ! must be false for 1D case !!
+
+   ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-)
+   LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h
+
+   ! advect (or not) open water. If not, retrieve it from advection of A
+   LOGICAL :: ll_ADVopw = .FALSE.  !
+   
+   ! 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
+   LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => true is better for 1D but false is better in 3D (for h and negative values) => pb in x-y?
+   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11 (worth than zalesak)
+
+   ! iterate on the limiter (only nonosc)
+   LOGICAL :: ll_limiter_it2 = .FALSE.
+   
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -39,10 +88,10 @@
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE ice_dyn_adv_umx( k_order, kt, pu_ice, pv_ice,  &
-      &                    pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ice_dyn_adv_umx  ***
@@ -54,5 +103,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the scheme (1-5 or 20)
+      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2)
       INTEGER                     , INTENT(in   ) ::   kt         ! time step
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
@@ -70,12 +119,14 @@
       !
       INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
-      INTEGER  ::   initad                  ! number of sub-timestep for the advection
-      REAL(wp) ::   zcfl , zusnit, zdt      !   -      -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box
+      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) ::   zcfl , zdt
+      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho
+      REAL(wp), DIMENSION(jpi,jpj) ::   zhvar
+      REAL(wp), DIMENSION(jpi,jpj) ::   zai_b, zai_a, z1_ai_b
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
       !
-      ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )
       !
       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
@@ -84,9 +135,9 @@
       IF( lk_mpp )   CALL mpp_max( zcfl )
 
-      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
-      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
-      ENDIF
-
-      zdt = rdt_ice / REAL(initad)
+      IF( zcfl > 0.5 ) THEN   ;   icycle = 2 
+      ELSE                    ;   icycle = 1 
+      ENDIF
+      
+      zdt = rdt_ice / REAL(icycle)
 
       ! --- transport --- !
@@ -109,34 +160,117 @@
       END DO
 
+      IF(.NOT. ALLOCATED(z1_ai))       ALLOCATE(z1_ai(jpi,jpj))
+      IF( ll_zeroup2 ) THEN
+         IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj))
+         IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj))
+      ENDIF
       !---------------!
       !== advection ==!
       !---------------!
-      DO jt = 1, initad
-         CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area 
+      DO jt = 1, icycle
+
+         IF( ll_ADVopw ) THEN
+            ll_dens=.FALSE.
+            zamsk = 1._wp
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area 
+         ELSE
+            zai_b(:,:) = SUM( pa_i(:,:,:), dim=3 )
+         ENDIF
+         
          DO jl = 1, jpl
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,jl) )        ! Salt content
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,jl) )        ! Age content
+            !
+            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_i(:,:,jl)
+            ELSEWHERE                         ;   z1_ai_b(:,:) = 0.
+            END WHERE
+            !
+            IF( ll_zeroup2 ) THEN
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), &
+                        &                              pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) )
+                     amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), &
+                        &                              pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) )
+                 END DO
+               END DO
+               CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.)
+            ENDIF
+            !
+            zamsk = 1._wp
+            ll_dens=.TRUE.
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ho, zva_ho ) ! Ice area
+            ll_dens=.FALSE.
+
+            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_i(:,:,jl)
+            ELSEWHERE                         ;   z1_ai(:,:) = 0.
+            END WHERE              
+
+            IF( ll_thickness ) THEN
+               zua_ho(:,:) = zudy(:,:)
+               zva_ho(:,:) = zvdx(:,:)
+            ENDIF
+            
+            zamsk = 0._wp ; zhvar(:,:) = pv_i (:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_i (:,:,jl) )              ! Ice volume
+            IF( ll_thickness )   pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
+            
+            zamsk = 0._wp ; zhvar(:,:) = pv_s (:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_s (:,:,jl) )              ! Snw volume
+            IF( ll_thickness )   pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
+
+            zamsk = 0._wp ; zhvar(:,:) = psv_i(:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), psv_i(:,:,jl) )              ! Salt content
+
+            zamsk = 0._wp ; zhvar(:,:) = poa_i(:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), poa_i(:,:,jl) )              ! Age content
+
             DO jk = 1, nlay_i
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content
-            END DO
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
+               zamsk = 0._wp ; zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai_b(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_i(:,:,jk,jl) )           ! Ice heat content
+            END DO
+
             DO jk = 1, nlay_s
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) )   ! Snow heat content
-            END DO
-            IF ( ln_pnd_H12 ) THEN
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume
-            ENDIF
-         END DO
-      END DO
-      !
-      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box )
+               zamsk = 0._wp ; zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai_b(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content
+            END DO
+            !
+            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_ai_b...)
+               !
+               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_ip(:,:,jl)
+               ELSEWHERE                          ;   z1_ai_b(:,:) = 0.
+               END WHERE
+               !
+               zamsk = 1._wp
+               ll_dens=.TRUE.
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ho, zva_ho ) ! mp fractio!n
+               ll_dens=.FALSE.
+
+               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_ip(:,:,jl)
+               ELSEWHERE                          ;   z1_ai(:,:) = 0.
+               END WHERE              
+
+               zamsk = 0._wp ; zhvar(:,:) = pv_ip(:,:,jl) * z1_ai_b(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_ip(:,:,jl) )              ! mp volume
+            ENDIF
+            !
+            !
+         END DO
+
+         IF( .NOT. ll_ADVopw ) THEN
+            zai_a(:,:) = SUM( pa_i(:,:,:), dim=3 )
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  pato_i(ji,jj) = pato_i(ji,jj) - ( zai_a(ji,jj) - zai_b(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( pato_i(:,:), 'T',  1. )
+         ENDIF
+         
+      END DO
       !
    END SUBROUTINE ice_dyn_adv_umx
 
    
-   SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc )
+   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -151,96 +285,639 @@
       !! ** Action : - pt  the after advective tracer
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme
-      INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration
-      REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptc            ! tracer content field
+      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(jpi,jpj), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt             ! tracer field
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
       !
       INTEGER  ::   ji, jj           ! dummy loop indices  
       REAL(wp) ::   ztra             ! local scalar
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfu_ho, zt_u, zt_ups
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfv_ups, zfv_ho, zt_v, ztrd
-      !!----------------------------------------------------------------------
-      !
-      !  upstream advection with initial mass fluxes & intermediate update
-      ! --------------------------------------------------------------------
-      DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj)
-            zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1)
-         END DO
-      END DO
+      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc 
+      !!----------------------------------------------------------------------
+      !
+      !  upstream (_ups) advection with initial mass fluxes
+      ! ---------------------------------------------------
+      IF( ll_clem )    zfu_ups=0.; zfv_ups=0.
+
+      IF( ll_gurvan .AND. pamsk==0. ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) &
+                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+            END DO
+         END DO
+         CALL lbc_lnk( pt, 'T', 1. )
+      ENDIF
+
       
-      DO jj = 2, jpjm1            ! total intermediate advective trends
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ztra = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
-               &       + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1)   ) * r1_e1e2t(ji,jj)
-            !
-            ztrd(ji,jj) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]  
-            zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)      ! guess after content field with monotonic scheme
-         END DO
-      END DO
-      CALL lbc_lnk( zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
-      
+      IF( .NOT. ll_upsxy ) THEN
+
+         ! fluxes
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               IF( ll_clem ) THEN
+                  zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
+               ELSE
+                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
+               ENDIF
+            END DO
+         END DO
+
+      ELSE
+         ! 1 if advection of A
+         ! z1_ai already defined IF advection of other tracers
+         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp
+         !
+         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+            ! flux in x-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( ll_clem ) THEN
+                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  ELSE
+                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  ENDIF
+               END DO
+            END DO
+            
+            ! first guess of tracer content from u-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
+                           &         * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) &
+                        &         * tmask(ji,jj,1)
+                  ENDIF
+!!                  IF( ji==26 .AND. jj==86) THEN
+!!                     WRITE(numout,*) '************************'
+!!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
+!!                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( zpt, 'T', 1. )
+            !
+            ! flux in y-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( ll_clem ) THEN
+                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1)
+                  ELSE
+                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
+                  ENDIF
+               END DO
+            END DO
+
+         !
+         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
+            ! flux in y-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( ll_clem ) THEN
+                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
+                  ELSE
+                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
+                  ENDIF
+               END DO
+            END DO
+
+            ! first guess of tracer content from v-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
+                        &                        + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
+                        &         * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
+                        &         * tmask(ji,jj,1)
+                  ENDIF
+!!                  IF( ji==26 .AND. jj==86) THEN
+!!                     WRITE(numout,*) '************************'
+!!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
+!!                  ENDIF
+                END DO
+            END DO
+            CALL lbc_lnk( zpt, 'T', 1. )
+            !
+            ! flux in x-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( ll_clem ) THEN
+                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj)
+                  ELSE
+                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
+                  ENDIF
+               END DO
+            END DO
+            !
+         ENDIF
+         
+      ENDIF
+
+      IF( ll_clem .AND. kn_limiter /= 1 )  &
+         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' )
+
+      IF( ll_zeroup2 ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( amaxu(ji,jj) == 0._wp )   zfu_ups(ji,jj) = 0._wp
+               IF( amaxv(ji,jj) == 0._wp )   zfv_ups(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+
+      ! guess after content field with upstream scheme
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
+               &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj)
+            IF( ll_clem ) THEN
+               IF( ll_gurvan ) THEN
+                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+               ELSE
+                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                                      + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
+            ELSE
+               zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+            ENDIF
+!!            IF( ji==26 .AND. jj==86) THEN
+!!               WRITE(numout,*) '**************************'
+!!               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj)
+!!            ENDIF
+         END DO
+      END DO
+      CALL lbc_lnk( zt_ups, 'T', 1. )
+
       ! High order (_ho) fluxes 
       ! -----------------------
-      SELECT CASE( k_order )
-      CASE ( 20 )                          ! centered second order
+      SELECT CASE( kn_umx )
+         !
+      CASE ( 20 )                          !== centered second order ==!
+         !
+         CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
+            &       zt_ups, zfu_ups, zfv_ups )
+         !
+      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
+         !
+         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
+            &        zt_ups, zfu_ups, zfv_ups )
+         !
+      END SELECT
+
+      IF( ll_thickness ) THEN
+         ! final trend with corrected fluxes
+         ! ------------------------------------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               IF( ll_gurvan ) THEN
+                  ztra       = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 
+               ELSE
+                  ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  & 
+                     &           + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) &
+                     &           + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj)
+               ENDIF
+               pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+
+               IF( pt(ji,jj) < -epsi20 ) THEN
+                  WRITE(numout,*) 'T<0 ',pt(ji,jj)
+               ENDIF
+               
+               IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 )   pt(ji,jj) = 0._wp
+               
+!!               IF( ji==26 .AND. jj==86) THEN
+!!                  WRITE(numout,*) 'zt high order',pt(ji,jj)
+!!               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk( pt, 'T',  1. )
+      ENDIF
+   
+      ! Rachid trick
+      ! ------------
+     IF( ll_clem ) THEN
+         IF( pamsk == 0. ) THEN
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
+                     zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj)
+                     zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj)
+                  ELSE
+                     zfu_ho (ji,jj) = 0._wp
+                     zfu_ups(ji,jj) = 0._wp
+                  ENDIF
+                  !
+                  IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
+                     zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj)
+                     zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj)
+                  ELSE
+                     zfv_ho (ji,jj) = 0._wp  
+                     zfv_ups(ji,jj) = 0._wp  
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDIF
+      ENDIF
+
+      IF( ll_zeroup5 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, fs_jpim1   ! vector opt.
+               zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                      - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1)
+               IF( zpt(ji,jj) < 0. ) THEN
+                  zfu_ho(ji,jj)   = zfu_ups(ji,jj)
+                  zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj)
+                  zfv_ho(ji,jj)   = zfv_ups(ji,jj)
+                  zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1)
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. )
+      ENDIF
+
+      ! output upstream trend of concentration and high order fluxes
+      ! ------------------------------------------------------------
+      IF( ll_dens ) THEN
+         ! high order u*a
          DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) )
-               zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) )
-            END DO
-         END DO
-         !
-      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme
-         CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
-         !
+            DO ji = 1, fs_jpim1
+               pua_ho (ji,jj) = zfu_ho (ji,jj)
+               pva_ho (ji,jj) = zfv_ho (ji,jj)
+            END DO
+         END DO
+         !!CALL lbc_lnk( pua_ho, 'U',  -1. ) ! clem: not needed I think
+         !!CALL lbc_lnk( pva_ho, 'V',  -1. )
+      ENDIF
+
+
+      IF( .NOT.ll_thickness ) THEN
+         ! final trend with corrected fluxes
+         ! ------------------------------------
          DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj)
-            END DO
-         END DO
+            DO ji = fs_2, fs_jpim1 
+               ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
+                  &     * r1_e1e2t(ji,jj) * pdt  
+
+               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN
+               !!   ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
+               !!      &     * r1_e1e2t(ji,jj) * pdt                 
+               !!ENDIF
+               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN
+               !!   WRITE(numout,*) 'Tc<0 ',ptc(ji,jj)+ztra
+               !!   ztra = 0._wp
+               !!ENDIF
+
+               ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1)
+                             
+!!               IF( ji==26 .AND. jj==86) THEN
+!!                  WRITE(numout,*) 'ztc high order',ptc(ji,jj)
+!!               ENDIF
+               
+            END DO
+         END DO
+         CALL lbc_lnk( ptc, 'T',  1. )
+      ENDIF
+      
+      !
+   END SUBROUTINE adv_umx
+
+   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
+      &             pt_ups, pfu_ups, pfv_ups )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE macho  ***
+      !!     
+      !! **  Purpose :   compute  
+      !!
+      !! **  Method  :   ... ???
+      !!                 TIM = transient interpolation Modeling 
+      !!
+      !! 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_limiter       ! limiter
+      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(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes 
+      !
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      LOGICAL  ::   ll_xy = .TRUE.
+      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --!
+         !
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj)
-            END DO
-         END DO
-         !
-      END SELECT
+            DO ji = 1, fs_jpim1
+               pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
+               pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
+            END DO
+         END DO
+         IF    ( kn_limiter == 1 ) THEN
+            IF( ll_clem ) THEN
+               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+            ELSE
+               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+            ENDIF
+         ELSEIF( kn_limiter == 2 ) THEN
+            CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
+            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
+         ELSEIF( kn_limiter == 3 ) THEN
+            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
+            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
+         ENDIF
+         !
+      ELSE                    !-- alternate directions --!
+         !
+         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp
+         !
+         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+            !
+            ! flux in x-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
+               END DO
+            END DO
+            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
+            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
+
+            ! first guess of tracer content from u-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1)
+                  ELSE                     
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+
+            ! flux in y-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
+               END DO
+            END DO
+            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
+            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
+
+         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
+            !
+            ! flux in y-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
+               END DO
+            END DO
+            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
+            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
+            !
+            ! first guess of tracer content from v-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1)
+                  ELSE
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+            !
+            ! flux in x-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
+               END DO
+            END DO
+            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
+            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
+
+         ENDIF
+         IF( ll_clem ) THEN
+            IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         ELSE
+            IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         ENDIF
          
-      ! antidiffusive flux : high order minus low order
-      ! --------------------------------------------------
-      DO jj = 2, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj)
-         END DO
-      END DO
-      DO jj = 1, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj)
-         END DO
-      END DO
-      
-      ! monotonicity algorithm
-      ! -------------------------
-      CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt )
-      
-      ! final trend with corrected fluxes
-      ! ------------------------------------
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.  
-            ztra       = ztrd(ji,jj)  - (  zfu_ho(ji,jj) - zfu_ho(ji-1,jj  )   &
-               &                         + zfv_ho(ji,jj) - zfv_ho(ji  ,jj-1) ) * r1_e1e2t(ji,jj)  
-            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      CALL lbc_lnk( ptc, 'T',  1. )
-      !
-   END SUBROUTINE adv_umx
-
-
-   SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
+      ENDIF
+   
+   END SUBROUTINE cen2
+
+   
+   SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, &
+      &              pt_ups, pfu_ups, pfv_ups )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE macho  ***
+      !!     
+      !! **  Purpose :   compute  
+      !!
+      !! **  Method  :   ... ???
+      !!                 TIM = transient interpolation Modeling 
+      !!
+      !! 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_limiter       ! limiter
+      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(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes 
+      !
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      REAL(wp) ::   ztra
+      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho
+      !!----------------------------------------------------------------------
+      !
+      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+         !
+         !                                                        !--  ultimate interpolation of pt at u-point  --!
+         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
+         !                                                        !--  limiter in x --!
+         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
+         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
+         !                                                        !--  advective form update in zzt  --!
+
+         IF( ll_1stguess_clem ) THEN
+
+            ! first guess of tracer content from u-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+
+         ELSE
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_gurvan ) THEN
+                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj)
+                  ELSE
+                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
+                  ENDIF
+                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+         ENDIF
+         !
+         !                                                        !--  ultimate interpolation of pt at v-point  --!
+         IF( ll_hoxy ) THEN
+            CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
+         ELSE
+            CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
+         ENDIF
+         !                                                        !--  limiter in y --!
+         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
+         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
+         !         
+         !
+      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
+         !
+         !                                                        !--  ultimate interpolation of pt at v-point  --!
+         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
+         !                                                        !--  limiter in y --!
+         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
+         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
+         !                                                        !--  advective form update in zzt  --!
+         IF( ll_1stguess_clem ) THEN
+            
+            ! first guess of tracer content from v-flux 
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
+                           &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
+                        &         * tmask(ji,jj,1)
+                  ENDIF
+                END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+            
+         ELSE
+            
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  IF( ll_gurvan ) THEN
+                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj)
+                  ELSE
+                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
+                  ENDIF
+                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+         ENDIF
+         !
+         !                                                        !--  ultimate interpolation of pt at u-point  --!
+         IF( ll_hoxy ) THEN
+            CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
+         ELSE
+            CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
+         ENDIF
+         !                                                        !--  limiter in x --!
+         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
+         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
+         !
+         !
+      ENDIF
+
+     
+      IF( kn_limiter == 1 ) THEN
+         IF( .NOT. ll_limiter_it2 ) THEN
+            IF( ll_clem ) THEN
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+            ELSE
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+            ENDIF
+         ELSE
+            zzfu_ho(:,:) = pfu_ho(:,:)
+            zzfv_ho(:,:) = pfv_ho(:,:)
+            ! 1st iteration of nonosc (limit the flux with the upstream solution)
+            IF( ll_clem ) THEN
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
+            ELSE
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
+            ENDIF
+            ! guess after content field with high order
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)
+                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+            ! 2nd iteration of nonosc (limit the flux with the limited high order solution)
+            IF( ll_clem ) THEN
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
+            ELSE
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+   END SUBROUTINE macho
+
+
+   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -253,78 +930,15 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
-      INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration
-      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points 
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zc_box    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zzt
-      !!----------------------------------------------------------------------
-      !
-      IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN         !==  odd ice time step:  adv_x then adv_y  ==!
-         !
-         !                                                           !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( k_order, pdt, ptc, puc, pt_u )
-         !
-         !                                                           !--  advective form update in zzt  --!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
-                  &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj)
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
-            END DO
-         END DO
-         CALL lbc_lnk( zzt, 'T', 1. )
-         !
-         !                                                           !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v )
-         !
-      ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==!
-         !
-         !                                                           !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v )
-         !
-         !                                                           !--  advective form update in zzt  --!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
-                  &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj)
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
-            END DO
-         END DO
-         CALL lbc_lnk( zzt, 'T', 1. )
-         !
-         !                                                           !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( k_order, pdt, zzt, puc, pt_u )
-         !      
-      ENDIF      
-      !
-   END SUBROUTINE macho
-
-
-   SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE ultimate_x  ***
-      !!     
-      !! **  Purpose :   compute  
-      !!
-      !! **  Method  :   ... ???
-      !!                 TIM = transient interpolation Modeling 
-      !!
-      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
-      !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
+      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
       REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity * A component
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point 
-      !
-      INTEGER  ::   ji, jj       ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho    ! high order flux 
+      !
+      INTEGER  ::   ji, jj             ! dummy loop indices
       REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztu1, ztu2, ztu3, ztu4
       !!----------------------------------------------------------------------
       !
@@ -354,12 +968,12 @@
       !
       !
-      SELECT CASE (k_order )
+      SELECT CASE (kn_umx )
       !
       CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
          !        
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   &
-                  &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
+               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
+                  &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
             END DO
          END DO
@@ -367,7 +981,7 @@
       CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
                   &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
@@ -377,7 +991,7 @@
       CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
 !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
@@ -391,7 +1005,7 @@
       CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
 !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
@@ -405,7 +1019,7 @@
       CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
 !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
@@ -421,9 +1035,30 @@
          !
       END SELECT
+      !                                                     !-- High order flux in i-direction  --!
+      IF( ll_neg ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               IF( pt_u(ji,jj) < 0._wp ) THEN
+                  pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
+                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
+               ENDIF
+            END DO
+         END DO
+      ENDIF
+
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            IF( ll_clem ) THEN
+               pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj)
+            ELSE
+               pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
+            ENDIF
+         END DO
+      END DO
       !
    END SUBROUTINE ultimate_x
    
  
-   SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v )
+   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_y  ***
@@ -436,13 +1071,15 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
+      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
       REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity*A component
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfv_ho    ! high order flux 
       !
       INTEGER  ::   ji, jj       ! dummy loop indices
       REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztv1, ztv2, ztv3, ztv4
       !!----------------------------------------------------------------------
       !
@@ -474,11 +1111,11 @@
       !
       !
-      SELECT CASE (k_order )
+      SELECT CASE (kn_umx )
       !
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  &
-                  &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
+            DO ji = 1, fs_jpim1
+               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
+                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
             END DO
          END DO
@@ -486,6 +1123,6 @@
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+            DO ji = 1, fs_jpim1
+               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
                   &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
@@ -496,6 +1133,6 @@
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+            DO ji = 1, fs_jpim1
+               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
 !!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
@@ -509,6 +1146,6 @@
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+            DO ji = 1, fs_jpim1
+               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
 !!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
@@ -522,6 +1159,6 @@
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+            DO ji = 1, fs_jpim1
+               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
 !!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
@@ -537,105 +1174,576 @@
          !
       END SELECT
+      !                                                     !-- High order flux in j-direction  --!
+      IF( ll_neg ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               IF( pt_v(ji,jj) < 0._wp ) THEN
+                  pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
+                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
+               ENDIF
+            END DO
+         END DO
+      ENDIF
+
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            IF( ll_clem ) THEN
+               pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj)
+            ELSE
+               pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
+            ENDIF
+         END DO
+      END DO
       !
    END SUBROUTINE ultimate_y
-   
-  
-   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt )
+     
+
+   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
       !!     
-      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
+       !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
       !!       scheme and the before field by a nonoscillatory algorithm 
       !!
       !! **  Method  :   ... ???
-      !!       warning : pbef and paft must be masked, but the boundaries
+      !!       warning : pt and pt_low must be masked, but the boundaries
       !!       conditions on the fluxes are not necessary zalezak (1979)
       !!       drange (1995) multi-dimensional forward-in-time and upstream-
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      REAL(wp)                     , INTENT(in   ) ::   pdt          ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pbef, paft   ! before & after field
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions
+      REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
+      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      INTEGER  ::   ikm1      ! local integer
-      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars
-      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv
-      !!----------------------------------------------------------------------
-      !
+      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt
+      !!----------------------------------------------------------------------
       zbig = 1.e+40_wp
-      zsml = 1.e-15_wp
-
-      ! test on divergence
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.  
-            zdiv(ji,jj) =  - (  paa(ji,jj) - paa(ji-1,jj  )   &
-               &              + pbb(ji,jj) - pbb(ji  ,jj-1) )  
-         END DO
-      END DO
-      CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
-
-      ! Determine ice masks for before and after tracers 
-      WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp )   ;   zmsk(:,:) = 0._wp
-      ELSEWHERE                                                                       ;   zmsk(:,:) = 1._wp * tmask(:,:,1)
-      END WHERE
-
+      zsml = epsi20
+
+      IF( ll_zeroup2 ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( amaxu(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp
+               IF( amaxv(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+      
+      IF( ll_zeroup4 ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( pfu_low(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp
+               IF( pfv_low(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+
+
+      IF( ll_zeroup1 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               IF( ll_gurvan ) THEN
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               ELSE
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
+               IF( zzt(ji,jj) < 0._wp ) THEN
+                  pfu_ho(ji,jj)   = pfu_low(ji,jj)
+                  pfv_ho(ji,jj)   = pfv_low(ji,jj)
+                  WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj)
+               ENDIF
+!!               IF( ji==26 .AND. jj==86) THEN
+!!                  WRITE(numout,*) 'zzt high order',zzt(ji,jj)
+!!                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+!!                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+!!                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
+!!                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
+!!               ENDIF
+               IF( ll_gurvan ) THEN
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               ELSE
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
+               IF( zzt(ji,jj) < 0._wp ) THEN
+                  pfu_ho(ji-1,jj) = pfu_low(ji-1,jj)
+                  pfv_ho(ji,jj-1) = pfv_low(ji,jj-1)
+                  WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj)
+               ENDIF
+               IF( ll_gurvan ) THEN
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               ELSE
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
+               IF( zzt(ji,jj) < 0._wp ) THEN
+                  WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj)
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )
+      ENDIF
+
+      
+      ! antidiffusive flux : high order minus low order
+      ! --------------------------------------------------
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj)
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj)
+          END DO
+      END DO
+
+      ! extreme case where pfu_ho has to be zero
+      ! ----------------------------------------
+      !                                    pfu_ho
+      !                           *         --->
+      !                        |      |  *   |        | 
+      !                        |      |      |    *   |    
+      !                        |      |      |        |    *
+      !            t_low :       i-1     i       i+1       i+2   
+      IF( ll_prelimiter_zalesak ) THEN
+         
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1 
+               zti_low(ji,jj)= pt_low(ji+1,jj  )
+               ztj_low(ji,jj)= pt_low(ji  ,jj+1)
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
+
+
+         !! this does not work
+         !!            DO jj = 2, jpjm1
+         !!               DO ji = fs_2, fs_jpim1
+         !!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj  ) - pt_low (ji  ,jj) ) .AND.     &
+         !!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj+1) - pt_low (ji  ,jj) )           &
+         !!               &    ) THEN
+         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji  ,jj) ) .AND.   &
+         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji  ,jj) )         &
+         !!               &       ) THEN
+         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0.
+         !!                     ENDIF
+         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj  ) ) .AND.  &
+         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji  ,jj-1) )        &
+         !!               &       )  THEN
+         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0.
+         !!                     ENDIF
+         !!                  ENDIF
+         !!                END DO
+         !!            END DO            
+
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND.  &
+                  & pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN
+                  !
+                  IF(  pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND.  &
+                     & pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0.
+                  !
+                  IF(  pfu_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji-1,jj) ) <= 0 .AND.  &
+                     & pfv_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji,jj-1) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0.
+                  !
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+
+      ELSEIF( ll_prelimiter_devore ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1 
+               zti_low(ji,jj)= pt_low(ji+1,jj  )
+               ztj_low(ji,jj)= pt_low(ji  ,jj+1)
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
+
+         z1_dt = 1._wp / pdt
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) )
+               pfu_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , &
+                  &                          zsign * ( pt_low (ji  ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji  ,jj) * z1_dt , &
+                  &                          zsign * ( zti_low(ji+1,jj) - zti_low(ji  ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) )
+
+               zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) )
+               pfv_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , &
+                  &                          zsign * ( pt_low (ji,jj  ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj  ) * z1_dt , &
+                  &                          zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) * e1e2t(ji,jj+1) * z1_dt ) )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+            
+      ENDIF
+         
+      
       ! Search local extrema
       ! --------------------
-      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
-!      zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   &
-!         &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  )
-!      zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   &
-!         &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  )
-      zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   &
-         &             paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  )
-      zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   &
-         &             paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  )
-
+      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF    ( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN
+               zbup(ji,jj) = -zbig
+               zbdo(ji,jj) =  zbig
+            ELSEIF( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) > 0._wp ) THEN
+               zbup(ji,jj) = pt_low(ji,jj)
+               zbdo(ji,jj) = pt_low(ji,jj)
+            ELSEIF( pt(ji,jj) > 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN
+               zbup(ji,jj) = pt(ji,jj)
+               zbdo(ji,jj) = pt(ji,jj)
+            ELSE
+               zbup(ji,jj) = MAX( pt(ji,jj) , pt_low(ji,jj) )
+               zbdo(ji,jj) = MIN( pt(ji,jj) , pt_low(ji,jj) )
+            ENDIF
+         END DO
+      END DO
+
+ 
       z1_dt = 1._wp / pdt
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
             !
-            zup  = MAX(   zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ),   &        ! search max/min in neighbourhood
-               &                       zbup(ji  ,jj-1), zbup(ji  ,jj+1)    )
-            zdo  = MIN(   zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ),   &
-               &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    )
-               !
-            zpos = MAX( 0., paa(ji-1,jj  ) ) - MIN( 0., paa(ji  ,jj  ) )   &        ! positive/negative  part of the flux
-               & + MAX( 0., pbb(ji  ,jj-1) ) - MIN( 0., pbb(ji  ,jj  ) )
-            zneg = MAX( 0., paa(ji  ,jj  ) ) - MIN( 0., paa(ji-1,jj  ) )   &
-               & + MAX( 0., pbb(ji  ,jj  ) ) - MIN( 0., pbb(ji  ,jj-1) )
-               !
-            zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms
-            zbetup(ji,jj) = ( zup         - paft(ji,jj) ) / ( zpos + zsml ) * zbt
-            zbetdo(ji,jj) = ( paft(ji,jj) - zdo         ) / ( zneg + zsml ) * zbt
+            IF( .NOT. ll_9points ) THEN
+            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1) )  ! search max/min in neighbourhood
+            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) )
+            !
+            ELSE
+            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1), &  ! search max/min in neighbourhood
+               &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  )
+            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), &
+               &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  )
+            ENDIF
+            !
+            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux
+               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )
+            zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) &
+               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) )
+            !
+            IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN
+               zneg2 = (   pt(ji,jj) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
+               zpos2 = ( - pt(ji,jj) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
+            ELSE
+               zneg2 = 0. ; zpos2 = 0.
+            ENDIF
+            !
+            !                                  ! up & down beta terms
+!            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
+!            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
+
+            IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt
+            ELSE                         ; zbetup(ji,jj) = 0. ! zbig
+            ENDIF
+            !
+            IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt
+            ELSE                         ; zbetdo(ji,jj) = 0. ! zbig
+            ENDIF
+            !
+            ! if all the points are outside ice cover
+            IF( zup == -zbig )   zbetup(ji,jj) = 0. ! zbig
+            IF( zdo ==  zbig )   zbetdo(ji,jj) = 0. ! zbig            
+            !
+
+!!            IF( ji==26 .AND. jj==86) THEN
+!               WRITE(numout,*) '-----------------'
+!               WRITE(numout,*) 'zpos',zpos,zpos2
+!               WRITE(numout,*) 'zneg',zneg,zneg2
+!               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj)))
+!               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj)))
+!               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj)))
+!               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1)))
+!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
+!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
+!               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
+!               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
+!               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt
+!               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt
+!               
+!               WRITE(numout,*) 'pt',pt(ji,jj)
+!               WRITE(numout,*) 'ptim1',pt(ji-1,jj)
+!               WRITE(numout,*) 'ptjm1',pt(ji,jj-1)
+!               WRITE(numout,*) 'pt_low',pt_low(ji,jj)
+!               WRITE(numout,*) 'zbetup',zbetup(ji,jj)
+!               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj)
+!               WRITE(numout,*) 'zup',zup
+!               WRITE(numout,*) 'zdo',zdo
+!            ENDIF
+            !
          END DO
       END DO
       CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
 
-      ! monotonic flux in the i & j direction (paa & pbb)
-      ! -------------------------------------
-      DO jj = 2, jpjm1
+      
+      ! monotonic flux in the y direction
+      ! ---------------------------------
+      DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector opt.
             zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
             zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
-            zcu = 0.5  + SIGN( 0.5 , paa(ji,jj) )
-            !
-            paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu )
-         END DO
-      END DO
-      !
+            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
+            !
+            zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
+            
+            pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj)
+
+!!            IF( ji==26 .AND. jj==86) THEN
+!!               WRITE(numout,*) 'coefU',zcoef
+!!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+!!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
+!!            ENDIF
+
+         END DO
+      END DO
+
       DO jj = 1, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO ji = 1, fs_jpim1   ! vector opt.
             zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
             zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
-            zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj) )
-            !
-            pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv )
-         END DO
-      END DO
+            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
+            !
+            zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
+            
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj)
+
+!!            IF( ji==26 .AND. jj==86) THEN
+!!               WRITE(numout,*) 'coefV',zcoef
+!!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+!!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
+!!            ENDIF
+         END DO
+      END DO
+
+      ! clem test
+      DO jj = 2, jpjm1
+         DO ji = 2, fs_jpim1   ! vector opt.
+            IF( ll_gurvan ) THEN
+               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+            ELSE
+               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+            ENDIF
+            IF( zzt(ji,jj) < -epsi20 ) THEN
+               WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj)
+            ENDIF
+         END DO
+      END DO
+      
+      !
       !
    END SUBROUTINE nonosc_2d
+
+   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE limiter_x  ***
+      !!     
+      !! **  Purpose :   compute flux limiter 
+      !!----------------------------------------------------------------------
+      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfu_ho       ! high order flux
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux
+      !
+      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      REAL(wp), DIMENSION (jpi,jpj) ::   zslpx       ! tracer slopes 
+      !!----------------------------------------------------------------------
+      !
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1)
+         END DO
+      END DO
+      CALL lbc_lnk( zslpx, 'U', -1.)   ! lateral boundary cond.
+      
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
+
+            Rjm = zslpx(ji-1,jj)
+            Rj  = zslpx(ji  ,jj)
+            Rjp = zslpx(ji+1,jj)
+
+            IF( PRESENT(pfu_ups) ) THEN
+
+               IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
+               ELSE                        ;   Rr = Rjp
+               ENDIF
+
+               zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj)     
+               IF( Rj > 0. ) THEN
+                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)),  &
+                     &        MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
+               ELSE
+                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)),  &
+                     &        MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
+               ENDIF
+               pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter
+               
+            ELSE
+               IF( Rj /= 0. ) THEN
+                  IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
+                  ELSE                        ;   Cr = Rjp / Rj
+                  ENDIF
+               ELSE
+                  Cr = 0.
+                  !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
+                  !ELSE                        ;   Cr = Rjp * 1.e20
+                  !ENDIF
+               ENDIF
+               
+               ! -- superbee --
+               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
+               ! -- van albada 2 --
+               !!zpsi = 2.*Cr / (Cr*Cr+1.)
+
+               ! -- sweby (with beta=1) --
+               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
+               ! -- van Leer --
+               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
+               ! -- ospre --
+               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
+               ! -- koren --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
+               ! -- charm --
+               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
+               !ELSE                 ;   zpsi = 0.
+               !ENDIF
+               ! -- van albada 1 --
+               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
+               ! -- smart --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
+               ! -- umist --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
+
+               ! high order flux corrected by the limiter
+               pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
+
+            ENDIF
+         END DO
+      END DO
+      CALL lbc_lnk( pfu_ho, 'U', -1.)   ! lateral boundary cond.
+      !
+   END SUBROUTINE limiter_x
+
+   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE limiter_y  ***
+      !!     
+      !! **  Purpose :   compute flux limiter 
+      !!----------------------------------------------------------------------
+      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfv_ho       ! high order flux
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux
+      !
+      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      REAL(wp), DIMENSION (jpi,jpj) ::   zslpy       ! tracer slopes 
+      !!----------------------------------------------------------------------
+      !
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1)
+         END DO
+      END DO
+      CALL lbc_lnk( zslpy, 'V', -1.)   ! lateral boundary cond.
+      
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
+
+            Rjm = zslpy(ji,jj-1)
+            Rj  = zslpy(ji,jj  )
+            Rjp = zslpy(ji,jj+1)
+
+            IF( PRESENT(pfv_ups) ) THEN
+
+               IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
+               ELSE                        ;   Rr = Rjp
+               ENDIF
+
+               zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj)     
+               IF( Rj > 0. ) THEN
+                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)),  &
+                     &        MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
+               ELSE
+                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)),  &
+                     &        MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
+               ENDIF
+               pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter
+
+            ELSE
+
+               IF( Rj /= 0. ) THEN
+                  IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
+                  ELSE                        ;   Cr = Rjp / Rj
+                  ENDIF
+               ELSE
+                  Cr = 0.
+                  !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
+                  !ELSE                        ;   Cr = Rjp * 1.e20
+                  !ENDIF
+               ENDIF
+
+               ! -- superbee --
+               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
+               ! -- van albada 2 --
+               !!zpsi = 2.*Cr / (Cr*Cr+1.)
+
+               ! -- sweby (with beta=1) --
+               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
+               ! -- van Leer --
+               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
+               ! -- ospre --
+               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
+               ! -- koren --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
+               ! -- charm --
+               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
+               !ELSE                 ;   zpsi = 0.
+               !ENDIF
+               ! -- van albada 1 --
+               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
+               ! -- smart --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
+               ! -- umist --
+               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
+
+               ! high order flux corrected by the limiter
+               pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
+
+            ENDIF
+         END DO
+      END DO
+      CALL lbc_lnk( pfv_ho, 'V', -1.)   ! lateral boundary cond.
+      !
+   END SUBROUTINE limiter_y
 
 #else
Index: NEMO/trunk/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10413)
@@ -39,15 +39,15 @@
 
    ! Variables shared among ridging subroutines
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   closing_net     ! net rate at which area is removed    (1/s)
-      !                                                 ! (ridging ice area - area of new ridges) / dt
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   opning         ! rate of opening due to divergence/shear
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   closing_gross  ! rate at which area removed, not counting area of new ridges
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmin    ! minimum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmax    ! maximum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hraft    ! thickness of rafted ice
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aridge   ! participating ice ridging
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   araft    ! participating ice rafting
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_net     ! net rate at which area is removed    (1/s)
+      !                                                               ! (ridging ice area - area of new ridges) / dt
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   opning          ! rate of opening due to divergence/shear
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_gross   ! rate at which area removed, not counting area of new ridges
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   apartf          ! participation function; fraction of ridging/closing associated w/ category n
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmin           ! minimum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmax           ! maximum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hraft           ! thickness of rafted ice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hi_hrdg         ! thickness of ridging ice / mean ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aridge          ! participating ice ridging
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   araft           ! participating ice rafting
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d
@@ -59,5 +59,4 @@
    LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
    REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
-   REAL(wp) ::   rn_crhg          ! determines changes in ice strength
    REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging            
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
@@ -79,5 +78,5 @@
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -193,5 +192,9 @@
             ! divergence given by the advection scheme
             !   (which may not be equal to divu as computed from the velocity field)
-            zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
+            IF    ( ln_adv_Pra ) THEN
+               zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
+            ELSEIF( ln_adv_UMx ) THEN
+               zdivu_adv(ji) = zdivu(ji)
+            ENDIF
             !
             IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
Index: NEMO/trunk/src/ICE/icedyn_rhg.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_rhg.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icedyn_rhg.F90	(revision 10413)
@@ -43,5 +43,5 @@
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
Index: NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 10413)
@@ -119,7 +119,9 @@
       REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
       REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
-      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
       REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
+      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
+      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
       !
       REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
@@ -137,4 +139,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param)
       REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
@@ -256,5 +259,9 @@
          END DO
       END DO
-            
+
+      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
+      IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile
+      ELSE                                               ;   zkt = 0._wp
+      ENDIF
       !
       !------------------------------------------------------------------------------!
@@ -317,4 +324,43 @@
       CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
       !
+      !                                  !== Landfast ice parameterization ==!
+      !
+      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ! ice thickness at U-V points
+               zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+               ! ice-bottom stress at U points
+               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
+               zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               ! ice-bottom stress at V points
+               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
+               zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               ! ice_bottom stress at T points
+               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+            END DO
+         END DO
+         CALL lbc_lnk( tau_icebfr(:,:), 'T', 1. )
+         !
+      ELSEIF( ln_landfast_home ) THEN          !-- Home made
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zTauU_ib(ji,jj) = tau_icebfr(ji,jj)
+               zTauV_ib(ji,jj) = tau_icebfr(ji,jj)
+            END DO
+         END DO
+         !
+      ELSE                                     !-- no landfast
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zTauU_ib(ji,jj) = 0._wp
+               zTauV_ib(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+      IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )
+
       !------------------------------------------------------------------------------!
       ! 3) Solution of the momentum equation, iterative procedure
@@ -380,7 +426,7 @@
                ENDIF
                
-               ! stress at T points
-               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
-               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
+               ! stress at T points (zkt/=0 if landfast)
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
              
             END DO
@@ -401,6 +447,6 @@
                zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
                
-               ! stress at F points
-               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
+               ! stress at F points (zkt/=0 if landfast)
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
 
             END DO
@@ -450,6 +496,6 @@
                   !
                   !                 !--- tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = - zTauV_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at V-points (energy conserving formulation)
@@ -462,5 +508,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
@@ -498,6 +544,6 @@
                   !
                   !                 !--- tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = - zTauU_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at U-points (energy conserving formulation)
@@ -510,5 +556,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
@@ -548,6 +594,6 @@
                   !
                   !                 !--- tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = - zTauU_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at U-points (energy conserving formulation)
@@ -560,5 +606,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
@@ -596,6 +642,6 @@
                   !
                   !                 !--- tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
-                  ztauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = - zTauV_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at v-points (energy conserving formulation)
@@ -608,5 +654,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
Index: NEMO/trunk/src/ICE/iceistate.F90
===================================================================
--- NEMO/trunk/src/ICE/iceistate.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/iceistate.F90	(revision 10413)
@@ -64,5 +64,5 @@
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -174,102 +174,110 @@
          ! then we check whether the distribution fullfills
          ! volume and area conservation, positivity and ice categories bounds
-         zh_i_ini(:,:,:) = 0._wp 
-         za_i_ini(:,:,:) = 0._wp
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
-
-                  ! find which category (jl0) the input ice thickness falls into
-                  jl0 = jpl
-                  DO jl = 1, jpl
-                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
-                        jl0 = jl
-                        CYCLE
-                     ENDIF
-                  END DO
+
+         IF( jpl == 1 ) THEN
+            !
+            zh_i_ini(:,:,1) = zht_i_ini(:,:)
+            za_i_ini(:,:,1) = zat_i_ini(:,:)            
+            !
+         ELSE
+            zh_i_ini(:,:,:) = 0._wp 
+            za_i_ini(:,:,:) = 0._wp
+            !
+            DO jj = 1, jpj
+               DO ji = 1, jpi
                   !
-                  itest(:) = 0
-                  i_fill   = jpl + 1                                            !------------------------------------
-                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
-                     !                                                          !------------------------------------
-                     i_fill = i_fill - 1
+                  IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
+
+                     ! find which category (jl0) the input ice thickness falls into
+                     jl0 = jpl
+                     DO jl = 1, jpl
+                        IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
+                           jl0 = jl
+                           CYCLE
+                        ENDIF
+                     END DO
                      !
-                     zh_i_ini(ji,jj,:) = 0._wp 
-                     za_i_ini(ji,jj,:) = 0._wp
                      itest(:) = 0
-                     !
-                     IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
-                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
-                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
-                     ELSE                         !-- case ice is thicker: fill categories >1
-                        ! thickness
-                        DO jl = 1, i_fill-1
-                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
-                        END DO
+                     i_fill   = jpl + 1                                            !------------------------------------
+                     DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
+                        !                                                          !------------------------------------
+                        i_fill = i_fill - 1
                         !
-                        ! concentration
-                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
-                        DO jl = 1, i_fill - 1
-                           IF( jl /= jl0 )THEN
-                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
-                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
+                        zh_i_ini(ji,jj,:) = 0._wp 
+                        za_i_ini(ji,jj,:) = 0._wp
+                        itest(:) = 0
+                        !
+                        IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
+                           zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
+                           za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
+                        ELSE                         !-- case ice is thicker: fill categories >1
+                           ! thickness
+                           DO jl = 1, i_fill-1
+                              zh_i_ini(ji,jj,jl) = hi_mean(jl)
+                           END DO
+                           !
+                           ! concentration
+                           za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
+                           DO jl = 1, i_fill - 1
+                              IF( jl /= jl0 )THEN
+                                 zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
+                                 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
+                              ENDIF
+                           END DO
+
+                           ! last category
+                           za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
+                           zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
+                           zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
+
+                           ! correction if concentration of upper cat is greater than lower cat
+                           !   (it should be a gaussian around jl0 but sometimes it is not)
+                           IF ( jl0 /= jpl ) THEN
+                              DO jl = jpl, jl0+1, -1
+                                 IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
+                                    zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
+                                    zh_i_ini(ji,jj,jl    ) = 0._wp
+                                    za_i_ini(ji,jj,jl    ) = 0._wp
+                                    za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
+                                       &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
+                                 END IF
+                              ENDDO
                            ENDIF
-                        END DO
-
-                        ! last category
-                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
-                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
-                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
-
-                        ! correction if concentration of upper cat is greater than lower cat
-                        !   (it should be a gaussian around jl0 but sometimes it is not)
-                        IF ( jl0 /= jpl ) THEN
-                           DO jl = jpl, jl0+1, -1
-                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
-                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
-                                 zh_i_ini(ji,jj,jl    ) = 0._wp
-                                 za_i_ini(ji,jj,jl    ) = 0._wp
-                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
-                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
-                              END IF
-                           ENDDO
+                           !
                         ENDIF
                         !
+                        ! Compatibility tests
+                        zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation
+                        IF ( zconv < epsi06 ) itest(1) = 1
+                        !
+                        zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation
+                           &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
+                        IF ( zconv < epsi06 ) itest(2) = 1
+                        !
+                        IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ?
+                        !
+                        itest(4) = 1
+                        DO jl = 1, i_fill
+                           IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations
+                        END DO
+                        !                                                          !----------------------------
+                     END DO                                                        ! end iteration on categories
+                     !                                                             !----------------------------
+                     IF( lwp .AND. SUM(itest) /= 4 ) THEN 
+                        WRITE(numout,*)
+                        WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
+                        WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure '
+                        WRITE(numout,*)
+                        WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
+                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
+                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
                      ENDIF
                      !
-                     ! Compatibility tests
-                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation
-                     IF ( zconv < epsi06 ) itest(1) = 1
-                     !
-                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation
-                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
-                     IF ( zconv < epsi06 ) itest(2) = 1
-                     !
-                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ?
-                     !
-                     itest(4) = 1
-                     DO jl = 1, i_fill
-                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations
-                     END DO
-                     !                                                          !----------------------------
-                  END DO                                                        ! end iteration on categories
-                  !                                                             !----------------------------
-                  IF( lwp .AND. SUM(itest) /= 4 ) THEN 
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
-                     WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure '
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
-                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
-                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
                   ENDIF
                   !
-               ENDIF
-               !
-            END DO   
-         END DO   
-
+               END DO
+            END DO
+         ENDIF
+         
          !---------------------------------------------------------------------
          ! 4) Fill in sea ice arrays
Index: NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- NEMO/trunk/src/ICE/icevar.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icevar.F90	(revision 10413)
@@ -557,4 +557,54 @@
       !!-------------------------------------------------------------------
       !
+      !
+      DO jl = 1, jpl       !==  loop over the categories  ==!
+         !
+         !----------------------------------------
+         ! zap ice energy and send it to the ocean
+         !----------------------------------------
+         DO jk = 1, nlay_i
+            DO jj = 1 , jpj
+               DO ji = 1 , jpi
+                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                     pe_i(ji,jj,jk,jl) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         DO jk = 1, nlay_s
+            DO jj = 1 , jpj
+               DO ji = 1 , jpi
+                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                     pe_s(ji,jj,jk,jl) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         !-----------------------------------------------------
+         ! zap ice and snow volume, add water and salt to ocean
+         !-----------------------------------------------------
+         DO jj = 1 , jpj
+            DO ji = 1 , jpi
+               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice
+                  pv_i   (ji,jj,jl) = 0._wp
+               ENDIF
+               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice
+                  pv_s   (ji,jj,jl) = 0._wp
+               ENDIF
+               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
+                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice
+                  psv_i  (ji,jj,jl) = 0._wp
+               ENDIF
+            END DO
+         END DO
+         !
+      END DO 
+      !
       WHERE( pato_i(:,:)   < 0._wp )   pato_i(:,:)   = 0._wp
       WHERE( poa_i (:,:,:) < 0._wp )   poa_i (:,:,:) = 0._wp
@@ -563,53 +613,4 @@
       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
       !                                                        but it does not change conservation, so keep it this way is ok
-      !
-      DO jl = 1, jpl       !==  loop over the categories  ==!
-         !
-         !----------------------------------------
-         ! zap ice energy and send it to the ocean
-         !----------------------------------------
-         DO jk = 1, nlay_i
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  IF( pe_i(ji,jj,jk,jl) < 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
-                     pe_i(ji,jj,jk,jl) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-         DO jk = 1, nlay_s
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  IF( pe_s(ji,jj,jk,jl) < 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
-                     pe_s(ji,jj,jk,jl) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-         !-----------------------------------------------------
-         ! zap ice and snow volume, add water and salt to ocean
-         !-----------------------------------------------------
-         DO jj = 1 , jpj
-            DO ji = 1 , jpi
-              IF( pv_i(ji,jj,jl) < 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice
-                  pv_i   (ji,jj,jl) = 0._wp
-               ENDIF
-               IF( pv_s(ji,jj,jl) < 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice
-                  pv_s   (ji,jj,jl) = 0._wp
-               ENDIF
-               IF( psv_i(ji,jj,jl) < 0._wp ) THEN
-                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice
-                  psv_i  (ji,jj,jl) = 0._wp
-               ENDIF
-            END DO
-         END DO
-         !
-      END DO 
       !
    END SUBROUTINE ice_var_zapneg
Index: NEMO/trunk/src/ICE/icewri.F90
===================================================================
--- NEMO/trunk/src/ICE/icewri.F90	(revision 10404)
+++ NEMO/trunk/src/ICE/icewri.F90	(revision 10413)
@@ -38,5 +38,5 @@
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
    !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -50,5 +50,5 @@
       INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::   z2da, z2db, zrho1, zrho2
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d !  2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast !  2D workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::   zmsk00, zmsk05, zmsk15, zmsksn ! O%, 5% and 15% concentration mask and snow mask
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zmsk00l, zmsksnl               ! cat masks
@@ -132,5 +132,6 @@
       IF( iom_use('vtau_ai'  ) )   CALL iom_put( "vtau_ai", vtau_ice * zmsk00     )   ! Wind stress term in force balance (y)
 
-      IF( iom_use('icevel') ) THEN 
+      IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN 
+        ! module of ice velocity
          DO jj = 2 , jpjm1
             DO ji = 2 , jpim1
@@ -141,5 +142,11 @@
          END DO
          CALL lbc_lnk( z2d, 'T', 1. )
-         IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d                   )   ! ice velocity module
+         IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d )
+
+        ! record presence of fast ice
+         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
+         ELSEWHERE                                                ; zfast(:,:) = 0._wp
+         END WHERE
+         IF( iom_use('fasticepres') )   CALL iom_put( "fasticepres" , zfast )
       ENDIF
 
Index: NEMO/trunk/src/NST/agrif_oce.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_oce.F90	(revision 10404)
+++ NEMO/trunk/src/NST/agrif_oce.F90	(revision 10413)
@@ -22,5 +22,5 @@
    !                                              !!* Namelist namagrif: AGRIF parameters
    LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !:
-   INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points)
+   INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 1  !: Sponge width (in number of parent grid points)
    REAL(wp), PUBLIC ::   rn_sponge_tra = 2800.     !: sponge coeff. for tracers
    REAL(wp), PUBLIC ::   rn_sponge_dyn = 2800.     !: sponge coeff. for dynamics
