Index: /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90	(revision 10266)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90	(revision 10267)
@@ -71,6 +71,8 @@
       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) ::   zhmax
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
       !!--------------------------------------------------------------------
       !
@@ -118,9 +120,35 @@
 
       CASE ( np_dynADV )           !==  pure advection ==!   (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 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)
+         ! --- monotonicity test from Lipscomb et al 2004 --- !
+         ! 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
+
+         ! 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
Index: /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90	(revision 10266)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90	(revision 10267)
@@ -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
@@ -43,6 +44,6 @@
 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 +55,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 +71,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) ::   zeps = 0.1_wp           ! shift in concentration to avoid division by 0
+      !                                     !    must be >= 0.01 and the best seems to be 0.1
+      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zfu, zfv
+      REAL(wp), DIMENSION(jpi,jpj) ::   z1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip 
       !!----------------------------------------------------------------------
       !
       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 +87,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
+      IF( zcfl > 0.5 ) THEN   ;   icycle = 2 
+      ELSE                    ;   icycle = 1 
       ENDIF
 
-      zdt = rdt_ice / REAL(initad)
+      zdt = rdt_ice / REAL(icycle)
 
       ! --- transport --- !
@@ -112,31 +115,76 @@
       !== 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
+
+         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 
+
          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
+            ! to avoid a problem with the limiter nonosc when A gets close to 0
+            pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1)
+            !
+            WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl)
+            ELSEWHERE                        ;   z1_a(:,:) = 0.   !!;   pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0. 
+            END WHERE
+            !
+            zamsk = 1._wp
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zfu, zfv ) ! Ice area
+            !!zfu = zudy ; zfv = zvdx
+            !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) )
+            zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl) )              ! Ice volume
+            zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) )              ! Snw volume
+            zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) )              ! Salt content
+            zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), 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 ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_i(:,:), 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
+               zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content
+            END DO
+            !
+            IF ( ln_pnd_H12 ) THEN   ! melt ponds
+               ! to avoid a problem with the limiter nonosc when A gets close to 0
+               pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1)
+               !
+               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_ap(:,:) = 1._wp / pa_ip(:,:,jl)
+               ELSEWHERE                         ;   z1_ap(:,:) = 0.
+               END WHERE
+               !
+               zamsk = 1._wp
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zfu, zfv ) ! mp fraction
+               zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl) )              ! mp volume
             ENDIF
-         END DO
-      END DO
-      !
-      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box )
+            !
+            ! to avoid a problem with the limiter nonosc when A gets close to 0
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1)
+                  pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps &
+                     &             + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
+                     &             ) * tmask(ji,jj,1)
+                  IF ( ln_pnd_H12 ) THEN   ! melt ponds
+                     pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps &
+                        &              + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
+                        &              ) * tmask(ji,jj,1)
+                  ENDIF
+              END DO
+            END DO
+            CALL lbc_lnk( pa_i(:,:,jl), 'T',  1. )
+            !
+         END DO
+
+      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, pfu, pfv )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -151,79 +199,86 @@
       !! ** 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(in   )           ::   pt             ! tracer field
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pfu, pfv       ! high order 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
-      
-      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)
+!!clem      REAL(wp) ::   zeps = 1.e-02        ! local scalar
+      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 
+      !!----------------------------------------------------------------------
+      !
+      ! add a constant value to avoid problems with zeros
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1)
+         END DO
+      END DO
+
+      !  upstream (_ups) advection with initial mass fluxes
+      ! ---------------------------------------------------
+      ! fluxes
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1
+            zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
+            zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
+         END DO
+      END DO
+      ! 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)
+            zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      CALL lbc_lnk( zt_ups, 'T', 1. )
       
       ! High order (_ho) fluxes 
       ! -----------------------
-      SELECT CASE( k_order )
-      CASE ( 20 )                          ! centered second order
-         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 jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj)
-               zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj)
-            END DO
-         END DO
+      SELECT CASE( kn_umx )
+         !
+      CASE ( 20 )                          !== centered second order ==!
+         !
+         CALL cen2( kn_limiter, jt, kt, pdt, zpt, 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, zpt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
+            &        zt_ups, zfu_ups, zfv_ups )
          !
       END SELECT
          
-      ! antidiffusive flux : high order minus low order
-      ! --------------------------------------------------
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj)
-            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 )
+      ! output high order fluxes
+      ! ------------------------
+      IF( PRESENT(pfu) ) THEN 
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               pfu(ji,jj) =  zfu_ho(ji,jj) !!clem - zeps * puc(ji,jj)
+               pfv(ji,jj) =  zfv_ho(ji,jj) !!clem - zeps * pvc(ji,jj)
+            END DO
+         END DO
+         !!CALL lbc_lnk( pfu, 'U',  -1. ) ! clem: not needed I think
+         !!CALL lbc_lnk( pfv, 'V',  -1. )
+      ENDIF
       
       ! 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)  
+         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)
+!!clem               &           + ( puc   (ji,jj) - puc   (ji-1,jj) + pvc   (ji,jj) - pvc   (ji,jj-1) ) * zeps  &  ! epsi * Div(ua)  or epsi * Div(u) 
+               &         ) * r1_e1e2t(ji,jj)  
             ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
          END DO
@@ -233,6 +288,200 @@
    END SUBROUTINE adv_umx
 
-
-   SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
+   SUBROUTINE cen2( 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. 
+      !!----------------------------------------------------------------------
+      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(in   ) ::   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 = 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
+            CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         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( 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.
+                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
+               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 * pv(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.
+                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
+               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 * pu(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( kn_limiter == 1 )   CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         
+      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(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes 
+      !
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
+      !!----------------------------------------------------------------------
+      !
+      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  --!
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               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
+               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( kn_umx, pdt, zzt, 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 )
+         !
+      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  --!
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               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
+               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( kn_umx, pdt, zzt, 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 )
+         !
+      ENDIF
+      IF( kn_limiter == 1 )   CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+      !
+   END SUBROUTINE macho
+
+
+   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -245,78 +494,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
       !!----------------------------------------------------------------------
       !
@@ -346,5 +532,5 @@
       !
       !
-      SELECT CASE (k_order )
+      SELECT CASE (kn_umx )
       !
       CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
@@ -352,6 +538,6 @@
          DO jj = 2, 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
@@ -361,5 +547,5 @@
          DO jj = 2, 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) ) ) 
@@ -371,5 +557,5 @@
          DO jj = 2, 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)
@@ -385,5 +571,5 @@
          DO jj = 2, 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)
@@ -399,5 +585,5 @@
          DO jj = 2, 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)
@@ -413,9 +599,15 @@
          !
       END SELECT
+      !                                                     !-- High order flux in i-direction  --!
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
+         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  ***
@@ -428,13 +620,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
       !!----------------------------------------------------------------------
       !
@@ -466,11 +660,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) ) )
+               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
@@ -479,5 +673,5 @@
          DO jj = 1, jpjm1
             DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               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) ) )
@@ -489,5 +683,5 @@
          DO jj = 1, jpjm1
             DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               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)
@@ -502,5 +696,5 @@
          DO jj = 1, jpjm1
             DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               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)
@@ -515,5 +709,5 @@
          DO jj = 1, jpjm1
             DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               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)
@@ -529,9 +723,15 @@
          !
       END SELECT
+      !                                                     !-- High order flux in j-direction  --!
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
+         END DO
+      END DO
       !
    END SUBROUTINE ultimate_y
-   
-  
-   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt )
+     
+
+   SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -541,20 +741,19 @@
       !!
       !! **  Method  :   ... ???
-      !!       warning : pbef and paft must be masked, but the boundaries
+      !!       warning : ptc and pt_ups 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   ) ::   pdt              ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_ups      ! before field & upstream guess of after field
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_ups, pfu_ups ! 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   ! local scalars
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo         !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv
+      !!----------------------------------------------------------------------
       zbig = 1.e+40_wp
       zsml = 1.e-15_wp
@@ -563,26 +762,32 @@
       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) )  
+            zdiv(ji,jj) =  - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) )  
          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
+      ! 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_ups(ji,jj)
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_ups(ji,jj)
+         END DO
+      END DO
 
       ! 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 ptc & pt_ups with large negative/positive value (-/+zbig) outside ice cover
+      DO jj = 1, jpj
+         DO ji = fs_2, fs_jpim1
+            IF( ptc(ji,jj) == 0._wp .AND. pt_ups(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN
+               zbup(ji,jj) = -zbig
+               zbdo(ji,jj) =  zbig
+            ELSE               
+               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_ups(ji,jj) )
+               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_ups(ji,jj) )
+            ENDIF
+         END DO
+      END DO
 
       z1_dt = 1._wp / pdt
@@ -590,44 +795,245 @@
          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
+            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) )
+            !
+            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &
+               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )  ! positive/negative part of the flux
+            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) )
+            !
+            !                                  ! up & down beta terms
+!!clem            zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
+!!clem            zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
+            IF( zpos >= epsi20 ) THEN
+               zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
+            ELSE
+               zbetup(ji,jj) = zbig
+            ENDIF
+            !
+            IF( zneg >= epsi20 ) THEN
+               zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
+            ELSE
+               zbetdo(ji,jj) = zbig
+            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) )
+            !
+            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj)
+         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 )
+            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
+            !
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj)
          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
