Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/ice.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/ice.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/ice.F90	(revision 10419)
@@ -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) **
@@ -325,10 +334,10 @@
    !! * Old values of global variables
    !!----------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b    !: snow and ice volumes/thickness
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b         !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                         !: snow heat content
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                         !: ice temperatures
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b              !: ice velocity
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                        !: ice concentration (total)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b, h_ip_b    !: snow and ice volumes/thickness
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b                 !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                                 !: snow heat content
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                                 !: ice temperatures
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b                      !: ice velocity
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                                !: ice concentration (total)
             
    !!----------------------------------------------------------------------
@@ -434,6 +443,6 @@
       ! * Old values of global variables
       ii = ii + 1
-      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   &
-         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,   &
+      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl), h_ip_b(jpi,jpj,jpl),  &
+         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,               &
          &      oa_i_b(jpi,jpj,jpl)                                                   , STAT=ierr(ii) )
 
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icectl.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icectl.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icectl.F90	(revision 10419)
@@ -197,7 +197,5 @@
 
       ! heat flux
-      zhfx  = glob_sum( 'icectl', ( qt_atm_oi - qt_oce_ai - diag_heat - diag_trp_ei - diag_trp_es   &
-      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
-         &              ) * e1e2t ) * zconv
+      zhfx  = glob_sum( 'icectl', ( qt_atm_oi - qt_oce_ai - diag_heat ) * e1e2t ) * zconv
 
       ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn.F90	(revision 10419)
@@ -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, zhip_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,126 @@
          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
+               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), &
+                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), &
+                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), &
+                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) )
+               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), &
+                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), &
+                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), &
+                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) )
+               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), &
+                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), &
+                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), &
+                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
+      !
+      !
+      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
+
+      CASE ( np_dynALL )           !==  all dynamical processes  ==!
+         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
+         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting 
+         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections
+
+      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
+         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology  
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
+         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting)
+
+      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( 'icedyn', 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( 'icedyn', 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, zhip_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( 'icedyn', 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, phip_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, phip_max   ! max ice thick from surrounding 9-pts
       !
       INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zh, zdv
+      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra
       !!-------------------------------------------------------------------
       !
@@ -156,14 +218,49 @@
          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_ip -- !
+                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
+                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN
+                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )
+                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN
+                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)
+                     ENDIF
                   ENDIF
                   !
+                  !                               ! -- check h_i -- !
+                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
+                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
+                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
+                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
+                  ENDIF
+                  !
+                  !                               ! -- check h_s -- !
+                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
+                  zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
+                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
+                     zfra = 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) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     !
+                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
+                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
+                  ENDIF           
+                  !
+                  !                               ! -- check snow load -- !
+                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
+                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
+                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
+                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+                  IF( zvs_excess > 0._wp ) THEN
+                     zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 )
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     !
+                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
+                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess
+                  ENDIF
+                  
                ENDIF
             END DO
@@ -215,7 +312,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 +330,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 +347,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 +363,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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv.F90	(revision 10419)
@@ -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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90	(revision 10419)
@@ -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,50 @@
    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
 
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: amaxu, amaxv
+   
+   ! 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 +86,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 +101,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,17 +117,22 @@
       !
       INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
-      INTEGER  ::   initad                  ! number of sub-timestep for the advection
-      INTEGER  ::   ipl                     ! third dimention of tracer array
-
-      REAL(wp) ::   zusnit, zdt
-      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zudy, zvdx, zcu_box, zcv_box
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zpato
+      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
+      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
+      REAL(wp) ::   zdt
+      REAL(wp), DIMENSION(1)       ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv
+      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx
+      REAL(wp), DIMENSION(jpi,jpj) ::   zati1, zati2
+
+
+
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zcu_box, zcv_box
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar
+
       !!----------------------------------------------------------------------
       !
       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) )
-      ALLOCATE( zpato(jpi,jpj,1) )
       !
       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
@@ -93,9 +145,9 @@
       CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
 
-      IF( zcflprv(1) > .5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
-      ELSE                         ;   initad = 1   ;   zusnit = 1.0_wp
+      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
+      ELSE                         ;   icycle = 1
       ENDIF
-
-      zdt = rdt_ice / REAL(initad)
+      
+      zdt = rdt_ice / REAL(icycle)
 
       ! --- transport --- !
@@ -118,52 +170,116 @@
       END DO
 
-      zpato (:,:,1) = pato_i(:,:)
-
+      IF( ll_zeroup2 ) THEN
+         IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj,jpl))
+         IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj,jpl))
+      ENDIF
       !---------------!
       !== advection ==!
       !---------------!
-      DO jt = 1, initad
+      DO jt = 1, icycle
+
+!!$         IF( ll_ADVopw ) THEN
+!!$            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 
+!!$            zamsk = 0._wp
+!!$         ELSE
+            zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
+!!$         ENDIF
          
-         l_full_nf_update = .FALSE.   ! false: disable full North fold update (performances)
-         CALL adv_umx( k_order, kt,   1, zdt, zudy, zvdx, zcu_box, zcv_box, zpato(:,:,1) )        ! Open water area 
-         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,:) )         ! Ice area
-         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,:) )         ! Ice  volume
-         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,:) )        ! Salt content
-         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,:) )        ! Age content
+         WHERE( pa_i(:,:,:) >= epsi20 )   ;   z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
+         ELSEWHERE                        ;   z1_ai(:,:,:) = 0.
+         END WHERE
+            !
+         WHERE( pa_ip(:,:,:) >= epsi20 )  ;   z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
+         ELSEWHERE                        ;   z1_aip(:,:,:) = 0.
+         END WHERE
+         !
+         IF( ll_zeroup2 ) THEN
+            DO jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     amaxu(ji,jj,jl)=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,jl)=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
+            END DO
+            CALL lbc_lnk_multi('icedyn_adv_umx', amaxu, 'T', 1., amaxv, 'T', 1.)
+         ENDIF
+         !
+         DO jl = 1, jpl
+            zua_ho(:,:,jl) = zudy(:,:)
+            zva_ho(:,:,jl) = zvdx(:,:)
+         END DO
+         
+         zamsk = 1._wp
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area
+         zamsk = 0._wp
+         !
+         IF( ll_thickness ) THEN
+            DO jl = 1, jpl
+               zua_ho(:,:,jl) = zudy(:,:)
+               zva_ho(:,:,jl) = zvdx(:,:)
+            END DO
+         ENDIF
+            !
+         zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_i )    ! Ice volume
+         IF( ll_thickness )   pv_i(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:)
+         !
+         zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_s )    ! Snw volume
+         IF( ll_thickness )   pv_s(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:)
+         !
+         zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, psv_i )    ! Salt content
+         !
+         zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, poa_i )    ! Age content
+         !
          DO jk = 1, nlay_i
-            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,:) )   ! Ice  heat content
-         END DO
-         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,:) )         ! Snow volume
+            zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content
+         END DO
+         !
          DO jk = 1, nlay_s
-            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,:) )   ! Snow heat content
-         END DO
+            zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content
+         END DO
+            !
          IF ( ln_pnd_H12 ) THEN
-            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,:) )     ! Melt pond fraction
-            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,:) )     ! Melt pond volume
+            !
+            DO jl = 1, jpl
+               zua_ho(:,:,jl) = zudy(:,:)
+               zva_ho(:,:,jl) = zvdx(:,:)
+            END DO
+            
+            zamsk = 1._wp
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction
+            zamsk = 0._wp
+            !
+            zhvar(:,:,:) = pv_ip(:,:,:) * z1_ai(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume
          ENDIF
-
-         l_full_nf_update = jt == initad   ! true: enable full North fold update (performances) when exiting the loop
-         CALL lbc_lnk( 'icedyn_adv_umx', zpato, 'T',  1. )
-         CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T',  1. )
-         CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T',  1. )
-         IF ( ln_pnd_H12 ) THEN
-            CALL lbc_lnk_multi( 'icedyn_adv_umx',  pa_i, 'T',  1., pv_i, 'T',  1., psv_i, 'T',  1., &
-                                                & poa_i, 'T',  1., pv_s, 'T',  1., pa_ip, 'T',  1., &
-                                                & pv_ip, 'T',  1. )
-         ELSE
-            CALL lbc_lnk_multi( 'icedyn_adv_umx',  pa_i, 'T',  1., pv_i, 'T',  1., psv_i, 'T',  1., &
-                                                & poa_i, 'T',  1., pv_s, 'T',  1. )
-         ENDIF
-         
+         !
+         !
+!!$         IF( .NOT. ll_ADVopw ) THEN
+            zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                  ! Open water area
+                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. )
+!!$         ENDIF
+         !
       END DO
       !
-      pato_i(:,:) = zpato (:,:,1)
-      !
-      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box, zpato )
-      !
    END SUBROUTINE ice_dyn_adv_umx
 
    
-   SUBROUTINE adv_umx( k_order, kt, ipl, 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  ***
@@ -178,112 +294,702 @@
       !! ** Action : - pt  the after advective tracer
       !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme
-      INTEGER                         , INTENT(in   ) ::   kt             ! number of iteration
-      INTEGER                         , INTENT(in   ) ::   ipl            ! third dimension of tracer array
-      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,ipl), 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(:,:  ), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
+      REAL(wp), DIMENSION(:,:  ), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   pt             ! tracer field
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   ptc            ! tracer content field
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
       !
       INTEGER  ::   ji, jj, jl       ! dummy loop indices  
-
       REAL(wp) ::   ztra             ! local scalar
-      REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfu_ups, zfu_ho, zt_u, zt_ups
-      REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfv_ups, zfv_ho, zt_v, ztrd
-
-      DO jl = 1, ipl
-      !!----------------------------------------------------------------------
-      !
-      !  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,jl) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj,jl)
-            zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1,jl)
+      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc 
+      !!----------------------------------------------------------------------
+      !
+      !  upstream (_ups) advection with initial mass fluxes
+      ! ---------------------------------------------------
+
+      IF( ll_gurvan .AND. pamsk==0. ) THEN
+         DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  pt(ji,jj,jl) = ( pt (ji,jj,jl) + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj)     &
+                     &                           + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) )   &
+                     &           * tmask(ji,jj,1)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. )
+      ENDIF
+
+      
+      IF( .NOT. ll_upsxy ) THEN
+
+         ! fluxes in both x-y directions
+         DO jl = 1, jpl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( ll_clem ) THEN
+                     zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
+                     zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
+                  ELSE
+                     zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl)
+                     zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+      ELSE
+         !
+         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 jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
+                     ELSE
+                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+            
+            ! first guess of tracer content from u-flux
+            DO jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     IF( ll_clem ) THEN
+                        IF( ll_gurvan ) THEN
+                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                        ELSE
+                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                              &            ) * tmask(ji,jj,1)
+                        ENDIF
+                     ELSE
+                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * 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
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+            !
+            ! flux in y-direction
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl)
+                     ELSE
+                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj+1,jl)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+         !
+         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
+            ! flux in y-direction
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
+                     ELSE
+                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+
+            ! first guess of tracer content from v-flux
+            DO jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     IF( ll_clem ) THEN
+                        IF( ll_gurvan ) THEN
+                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                        ELSE
+                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
+                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
+                              &            * tmask(ji,jj,1)
+                        ENDIF
+                     ELSE
+                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * 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
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
+            !
+            ! flux in x-direction
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl)
+                     ELSE
+                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * zpt(ji+1,jj,jl)
+                     ENDIF
+                  END DO
+               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 jl = 1, jpl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  IF( amaxu(ji,jj,jl) == 0._wp )   zfu_ups(ji,jj,jl) = 0._wp
+                  IF( amaxv(ji,jj,jl) == 0._wp )   zfv_ups(ji,jj,jl) = 0._wp
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      ! guess after content field with upstream scheme
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ztra          = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   &
+                  &                + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj)
+               IF( ll_clem ) THEN
+                  IF( ll_gurvan ) THEN
+                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
+                  ELSE
+                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) )   &
+                        &                                            +   pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) &
+                        &                                              * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+                  ENDIF
+               ELSE
+                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + 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
       END DO
-      
-         DO jj = 2, jpjm1            ! total intermediate advective trends
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ztra = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   &
-                  &       + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl)   ) * r1_e1e2t(ji,jj)
-            !
-               ztrd(ji,jj,jl) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]  
-               zt_ups (ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)   ! guess after content field with monotonic scheme
-            END DO
-         END DO
-      END DO
-      
+      CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )
+
       ! High order (_ho) fluxes 
       ! -----------------------
-      SELECT CASE( k_order )
-      CASE ( 20 )                          ! centered second order
-         DO jl = 1, ipl
+      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 jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  IF( ll_gurvan ) THEN
+                     ztra       = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 
+                  ELSE
+                     ztra       = ( - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  & 
+                        &           + ( pt(ji,jj,jl) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) &
+                        &           + ( pt(ji,jj,jl) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj)
+                  ENDIF
+                  pt(ji,jj,jl) = ( pt(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
+                  
+                  IF( pt(ji,jj,jl) < -epsi20 ) THEN
+                     WRITE(numout,*) 'T<0 ',pt(ji,jj,jl)
+                  ENDIF
+                  
+                  IF( pt(ji,jj,jl) < 0._wp .AND. pt(ji,jj,jl) >= -epsi20 )   pt(ji,jj,jl) = 0._wp
+                  
+                  !!               IF( ji==26 .AND. jj==86) THEN
+                  !!                  WRITE(numout,*) 'zt high order',pt(ji,jj)
+                  !!               ENDIF
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T',  1. )
+      ENDIF
+   
+      ! Rachid trick
+      ! ------------
+      IF( ll_clem ) THEN
+         IF( pamsk == 0. ) THEN
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
+                        zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
+                        zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
+                     ELSE
+                        zfu_ho (ji,jj,jl) = 0._wp
+                        zfu_ups(ji,jj,jl) = 0._wp
+                     ENDIF
+                     !
+                     IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
+                        zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
+                        zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
+                     ELSE
+                        zfv_ho (ji,jj,jl) = 0._wp  
+                        zfv_ups(ji,jj,jl) = 0._wp  
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+
+      IF( ll_zeroup5 ) THEN
+         DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = 2, fs_jpim1   ! vector opt.
+                  zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                            - ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1)
+                  IF( zpt(ji,jj,jl) < 0. ) THEN
+                     zfu_ho(ji  ,jj,jl) = zfu_ups(ji  ,jj,jl)
+                     zfu_ho(ji-1,jj,jl) = zfu_ups(ji-1,jj,jl)
+                     zfv_ho(ji  ,jj,jl) = zfv_ups(ji  ,jj,jl)
+                     zfv_ho(ji,jj-1,jl) = zfv_ups(ji,jj-1,jl)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. )
+      ENDIF
+
+      ! output high order fluxes u*a
+      ! ----------------------------
+      IF( PRESENT( pua_ho ) ) THEN
+         DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji+1,jj,jl) )
-                  zfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji,jj+1,jl) )
-               END DO
-            END DO
-         END DO
-         !
-      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme
-         CALL macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
-         !
-         DO jl = 1, ipl
+               DO ji = 1, fs_jpim1
+                  pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl)
+                  pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+
+      IF( .NOT.ll_thickness ) THEN
+         ! final trend with corrected fluxes
+         ! ------------------------------------
+         DO jl = 1, jpl
             DO jj = 2, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zfu_ho(ji,jj,jl) = puc(ji,jj) * zt_u(ji,jj,jl)
-               END DO
-            END DO
-         END DO
-         DO jl = 1, ipl
+               DO ji = fs_2, fs_jpim1 
+                  ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt  
+                  
+                  ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + 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
+         END DO
+         CALL lbc_lnk( 'icedyn_adv_umx', 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(:,:,:), INTENT(in   ) ::   pt               ! tracer fields
+      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step 
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes 
+      !
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
+      LOGICAL  ::   ll_xy = .TRUE.
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --!
+         !
+         DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zfv_ho(ji,jj,jl) = pvc(ji,jj) * zt_v(ji,jj,jl)
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
+               DO ji = 1, fs_jpim1
+                  IF( ll_clem ) THEN
+                     pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
+                     pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
+                  ELSE
+                     pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
+                     pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
+                  ENDIF
+               END DO
+            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( 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 jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
+                     ELSE
+                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
+                     ENDIF
+                  END DO
+               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 jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     IF( ll_clem ) THEN
+                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)        &
+                           &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
+                           &         * tmask(ji,jj,1)
+                     ELSE                     
+                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
+
+            ! flux in y-direction
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) )
+                     ELSE                     
+                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) )
+                     ENDIF
+                  END DO
+               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 jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
+                     ELSE                     
+                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
+                     ENDIF
+                  END DO
+               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 jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     IF( ll_clem ) THEN
+                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
+                           &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
+                           &         * tmask(ji,jj,1)
+                     ELSE
+                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
+            !
+            ! flux in x-direction
+            DO jl = 1, jpl
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1
+                     IF( ll_clem ) THEN
+                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) )
+                     ELSE
+                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) )
+                     ENDIF
+                  END DO
+               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 jl = 1, ipl
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zfu_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) - zfu_ups(ji,jj,jl)
-            END DO
-         END DO
-      END DO
-      DO jl = 1, ipl
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zfv_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) - zfv_ups(ji,jj,jl)
-            END DO
-         END DO
-      END DO
-
-      CALL lbc_lnk("icedyn_adv_umx",zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
-      
-      ! monotonicity algorithm
-      ! -------------------------
-      CALL nonosc_2d( ipl, ptc, zfu_ho, zfv_ho, zt_ups, pdt )
-      
-      ! final trend with corrected fluxes
-      ! ------------------------------------
-      DO jl = 1, ipl
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.  
-               ztra       = ztrd(ji,jj,jl)  - (  zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj  ,jl)   &
-                  &                         +    zfv_ho(ji,jj,jl) - zfv_ho(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj)  
-               ptc(ji,jj,jl) = ptc(ji,jj,jl) + pdt * ztra * tmask(ji,jj,1)
-            END DO
-         END DO
-      END DO
-      !
-   END SUBROUTINE adv_umx
-
-   SUBROUTINE macho( k_order, kt, ipl, 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(:,:,:), INTENT(in   ) ::   pt               ! tracer fields
+      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
+      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step 
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points 
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes 
+      !
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
+      REAL(wp) ::   ztra
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   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 jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     IF( ll_clem ) THEN
+                        IF( ll_gurvan ) THEN
+                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                        ELSE
+                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                              &         ) * tmask(ji,jj,1)
+                        ENDIF
+                     ELSE
+                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
+
+         ELSE
+
+            DO jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     IF( ll_gurvan ) THEN
+                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  &
+                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj)
+                     ELSE
+                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  &
+                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
+                     ENDIF
+                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', 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 jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     IF( ll_clem ) THEN
+                        IF( ll_gurvan ) THEN
+                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                        ELSE
+                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
+                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                              &         ) * tmask(ji,jj,1)
+                        ENDIF
+                     ELSE
+                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) &
+                           &         * tmask(ji,jj,1)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
+            
+         ELSE
+            
+            DO jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     IF( ll_gurvan ) THEN
+                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  &
+                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj)
+                     ELSE
+                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  &
+                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
+                     ENDIF
+                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', 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 jl = 1, jpl
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     ztra          = - ( zzfu_ho(ji,jj,jl) - zzfu_ho(ji-1,jj,jl) + zzfv_ho(ji,jj,jl) - zzfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj)
+                     zzt(ji,jj,jl) =   ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_adv_umx', 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  ***
@@ -296,96 +1002,26 @@
       !! 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
-      INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array
-      REAL(wp)                        , INTENT(in   ) ::   pdt        ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,ipl), 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,ipl), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points 
-      !
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      REAL(wp) ::   zc_box    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,ipl) :: 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, ipl, pdt, ptc, puc, pt_u )
-         !
-         !                                                           !--  advective form update in zzt  --!
-         DO jl = 1, ipl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zzt(ji,jj,jl) = ptc(ji,jj,jl) - pubox(ji,jj) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  &
-                   &            - ptc(ji,jj,jl) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj)
-                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
-         !
-         !                                                           !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( k_order, ipl, 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, ipl, pdt, ptc, pvc, pt_v )
-         !
-         !                                                           !--  advective form update in zzt  --!
-         DO jl = 1, ipl
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zzt(ji,jj,jl) = ptc(ji,jj,jl) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  &
-                     &                    - ptc  (ji,jj,jl) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj)
-                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
-         !
-         !                                                           !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( k_order, ipl, pdt, zzt, puc, pt_u )
-         !      
-      ENDIF      
-      !
-   END SUBROUTINE macho
-
-
-   SUBROUTINE ultimate_x( k_order, ipl, 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   ) ::   ipl        ! third dimension of tracer array
-      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,ipl), INTENT(in   ) ::   pt        ! tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u      ! tracer at u-point 
-      !
-      INTEGER  ::   ji, jj, jl       ! dummy loop indices
+      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(:,:  ), INTENT(in   ) ::   pu        ! ice i-velocity component
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc       ! ice i-velocity * A component
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u      ! tracer at u-point 
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho    ! high order flux 
+      !
+      INTEGER  ::   ji, jj, jl             ! dummy loop indices
       REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu3
-      REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztu2, ztu4
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4
       !!----------------------------------------------------------------------
       !
       !                                                     !--  Laplacian in i-direction  --!
-      DO jl = 1, ipl
+      DO jl = 1, jpl
          DO jj = 2, jpjm1         ! First derivative (gradient)
             DO ji = 1, fs_jpim1
-               ztu1(ji,jj) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+               ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
             END DO
             !                     ! Second derivative (Laplacian)
             DO ji = fs_2, fs_jpim1
-               ztu2(ji,jj,jl) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
+               ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
             END DO
          END DO
@@ -394,12 +1030,12 @@
       !
       !                                                     !--  BiLaplacian in i-direction  --!
-      DO jl = 1, ipl
+      DO jl = 1, jpl
          DO jj = 2, jpjm1         ! Third derivative
             DO ji = 1, fs_jpim1
-               ztu3(ji,jj) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+               ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
             END DO
             !                     ! Fourth derivative
             DO ji = fs_2, fs_jpim1
-               ztu4(ji,jj,jl) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
+               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
             END DO
          END DO
@@ -408,13 +1044,13 @@
       !
       !
-      SELECT CASE (k_order )
+      SELECT CASE (kn_umx )
       !
       CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
          !        
-         DO jl = 1, ipl
-            DO jj = 2, jpjm1
+         DO jl = 1, jpl
+            DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
-                     &                                       - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
+                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
                END DO
             END DO
@@ -423,10 +1059,10 @@
       CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
          !
-         DO jl = 1, ipl
-            DO jj = 2, jpjm1
+         DO jl = 1, jpl
+            DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
-                     &                                              -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
+                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
+                     &                                               -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
                END DO
             END DO
@@ -435,14 +1071,14 @@
       CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
          !
-         DO jl = 1, ipl
-            DO jj = 2, jpjm1
+         DO jl = 1, jpl
+            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)
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
-                     &                                              -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
-                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
-                     &                                              - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
+                     &                                               -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
+                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
+                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
                END DO
             END DO
@@ -451,14 +1087,14 @@
       CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
          !
-         DO jl = 1, ipl
-            DO jj = 2, jpjm1
+         DO jl = 1, jpl
+            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)
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
-                     &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
-                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
-                     &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
+                     &                                               -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
+                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
+                     &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
                END DO
             END DO
@@ -467,17 +1103,17 @@
       CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
          !
-         DO jl = 1, ipl
-            DO jj = 2, jpjm1
+         DO jl = 1, jpl
+            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)
+                  !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
                   zdx4 = zdx2 * zdx2
-                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (           (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       &
-                     &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   &
-                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *    (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       &
-                     &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   &
-                     &       + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       &
-                     &                                              - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (            (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       &
+                     &                                                     -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   &
+                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       &
+                     &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   &
+                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       &
+                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
                END DO
             END DO
@@ -485,9 +1121,34 @@
          !
       END SELECT
+      !                                                     !-- High order flux in i-direction  --!
+      IF( ll_neg ) THEN
+         DO jl = 1, jpl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN
+                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
+                        &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      DO jl = 1, jpl
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( ll_clem ) THEN
+                  pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
+               ELSE
+                  pfu_ho(ji,jj,jl) = puc(ji,jj,jl) * pt_u(ji,jj,jl)
+               ENDIF
+            END DO
+         END DO
+      END DO
       !
    END SUBROUTINE ultimate_x
    
  
-   SUBROUTINE ultimate_y( k_order, ipl, pdt, pt, pvc, pt_v )
+   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_y  ***
@@ -500,27 +1161,27 @@
       !! 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   ) ::   ipl       ! third dimension of tracer array
-      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,ipl), INTENT(in   ) ::   pt        ! tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_v      ! tracer at v-point 
+      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(:,:  ), INTENT(in   ) ::   pv        ! ice j-velocity component
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pvc       ! ice j-velocity*A component
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_v      ! tracer at v-point 
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfv_ho    ! high order flux 
       !
       INTEGER  ::   ji, jj, jl       ! dummy loop indices
       REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv3
-      REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztv2, ztv4
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4
       !!----------------------------------------------------------------------
       !
       !                                                     !--  Laplacian in j-direction  --!
-      DO jl = 1, ipl
+      DO jl = 1, jpl
          DO jj = 1, jpjm1         ! First derivative (gradient)
             DO ji = fs_2, fs_jpim1
-               ztv1(ji,jj) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
+               ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
             END DO
          END DO
          DO jj = 2, jpjm1         ! Second derivative (Laplacian)
             DO ji = fs_2, fs_jpim1
-               ztv2(ji,jj,jl) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
+               ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
             END DO
          END DO
@@ -529,13 +1190,13 @@
       !
       !                                                     !--  BiLaplacian in j-direction  --!
-      DO jl = 1, ipl
+      DO jl = 1, jpl
          DO jj = 1, jpjm1         ! First derivative
             DO ji = fs_2, fs_jpim1
-            ztv3(ji,jj) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
+               ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
             END DO
          END DO
          DO jj = 2, jpjm1         ! Second derivative
             DO ji = fs_2, fs_jpim1
-               ztv4(ji,jj,jl) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
+               ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
             END DO
          END DO
@@ -544,12 +1205,12 @@
       !
       !
-      SELECT CASE (k_order )
-      !
+      SELECT CASE (kn_umx )
+         !
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
-         DO jl = 1, ipl
+         DO jl = 1, jpl
             DO jj = 1, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
-                     &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
+               DO ji = 1, fs_jpim1
+                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
+                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
                END DO
             END DO
@@ -557,8 +1218,8 @@
          !
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
-         DO jl = 1, ipl
+         DO jl = 1, jpl
             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,jl) = 0.5_wp * vmask(ji,jj,1) * (     ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
                      &                                     - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -569,8 +1230,8 @@
          !
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
-         DO jl = 1, ipl
+         DO jl = 1, jpl
             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)
@@ -584,8 +1245,8 @@
          !
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
-         DO jl = 1, ipl
+         DO jl = 1, jpl
             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)
@@ -599,8 +1260,8 @@
          !
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
-         DO jl = 1, ipl
+         DO jl = 1, jpl
             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)
@@ -617,85 +1278,343 @@
          !
       END SELECT
+      !                                                     !-- High order flux in j-direction  --!
+      IF( ll_neg ) THEN
+         DO jl = 1, jpl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN
+                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
+                        &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      DO jl = 1, jpl
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( ll_clem ) THEN
+                  pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
+               ELSE
+                  pfv_ho(ji,jj,jl) = pvc(ji,jj,jl) * pt_v(ji,jj,jl)
+               ENDIF
+            END DO
+         END DO
+      END DO
       !
    END SUBROUTINE ultimate_y
-   
-  
-   SUBROUTINE nonosc_2d( ipl, 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
       !!----------------------------------------------------------------------
-      INTEGER                          , INTENT(in   ) ::   ipl          ! third dimension of tracer array
-      REAL(wp)                         , INTENT(in   ) ::   pdt          ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(in   ) ::   pbef, paft   ! before & after field
-      REAL(wp), DIMENSION (jpi,jpj,ipl), 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 (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a
+      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field
+      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux
+      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
       !
       INTEGER  ::   ji, jj, jl    ! 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) :: zbup, zbdo, zmsk
-      REAL(wp), DIMENSION(jpi,jpj,ipl) :: zbetup, zbetdo, 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    ) :: zbup, zbdo
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_low, ztj_low, zzt
+      !!----------------------------------------------------------------------
       zbig = 1.e+40_wp
-      zsml = 1.e-15_wp
-
-      ! test on divergence
-      DO jl = 1, ipl
+      zsml = epsi20
+      
+      IF( ll_zeroup2 ) THEN
+         DO jl = 1, jpl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  IF( amaxu(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp
+                  IF( amaxv(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      IF( ll_zeroup4 ) THEN
+         DO jl = 1, jpl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  IF( pfu_low(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp
+                  IF( pfv_low(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+
+      IF( ll_zeroup1 ) THEN
+         DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  IF( ll_gurvan ) THEN
+                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                  ELSE
+                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                        &         ) * tmask(ji,jj,1)
+                  ENDIF
+                  IF( zzt(ji,jj,jl) < 0._wp ) THEN
+                     pfu_ho(ji,jj,jl)   = pfu_low(ji,jj,jl)
+                     pfv_ho(ji,jj,jl)   = pfv_low(ji,jj,jl)
+                     WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)
+                  ENDIF
+                  !!               IF( ji==26 .AND. jj==86) THEN
+                  !!                  WRITE(numout,*) 'zzt high order',zzt(ji,jj)
+                  !!                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
+                  !!                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
+                  !!                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt
+                  !!                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt
+                  !!               ENDIF
+                  IF( ll_gurvan ) THEN
+                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                  ELSE
+                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                        &         ) * tmask(ji,jj,1)
+                  ENDIF
+                  IF( zzt(ji,jj,jl) < 0._wp ) THEN
+                     pfu_ho(ji-1,jj,jl) = pfu_low(ji-1,jj,jl)
+                     pfv_ho(ji,jj-1,jl) = pfv_low(ji,jj-1,jl)
+                     WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)
+                  ENDIF
+                  IF( ll_gurvan ) THEN
+                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                  ELSE
+                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                        &         ) * tmask(ji,jj,1)
+                  ENDIF
+                  IF( zzt(ji,jj,jl) < 0._wp ) THEN
+                     WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )
+      ENDIF
+
+      
+      ! antidiffusive flux : high order minus low order
+      ! --------------------------------------------------
+      DO jl = 1, jpl
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_low(ji,jj,jl)
+               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_low(ji,jj,jl)
+            END DO
+         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 jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1 
+                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl)
+                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', 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 jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  IF ( pfu_ho(ji,jj,jl) * ( pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) <= 0. .AND.  &
+                     & pfv_ho(ji,jj,jl) * ( pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) <= 0. ) THEN
+                     !
+                     IF(  pfu_ho(ji,jj,jl) * ( zti_low(ji+1,jj,jl) - zti_low(ji,jj,jl) ) <= 0 .AND.  &
+                        & pfv_ho(ji,jj,jl) * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0.
+                     !
+                     IF(  pfu_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji-1,jj,jl) ) <= 0 .AND.  &
+                        & pfv_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji,jj-1,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0.
+                     !
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+
+      ELSEIF( ll_prelimiter_devore ) THEN
+         DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1 
+                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl)
+                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. )
+
+         z1_dt = 1._wp / pdt
+         DO jl = 1, jpl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zsign = SIGN( 1., pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) )
+                  pfu_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj,jl)) , &
+                     &                          zsign * ( pt_low (ji  ,jj,jl) - pt_low (ji-1,jj,jl) ) * e1e2t(ji  ,jj) * z1_dt , &
+                     &                          zsign * ( zti_low(ji+1,jj,jl) - zti_low(ji  ,jj,jl) ) * e1e2t(ji+1,jj) * z1_dt ) )
+
+                  zsign = SIGN( 1., pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) )
+                  pfv_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj,jl)) , &
+                     &                          zsign * ( pt_low (ji,jj  ,jl) - pt_low (ji,jj-1,jl) ) * e1e2t(ji,jj  ) * z1_dt , &
+                     &                          zsign * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj  ,jl) ) * e1e2t(ji,jj+1) * z1_dt ) )
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+
+      ENDIF
+
+
+      ! Search local extrema
+      ! --------------------
+      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover
+      z1_dt = 1._wp / pdt
+      DO jl = 1, jpl
+         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN
+                  zbup(ji,jj) = -zbig
+                  zbdo(ji,jj) =  zbig
+               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) > 0._wp ) THEN
+                  zbup(ji,jj) = pt_low(ji,jj,jl)
+                  zbdo(ji,jj) = pt_low(ji,jj,jl)
+               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN
+                  zbup(ji,jj) = pt(ji,jj,jl)
+                  zbdo(ji,jj) = pt(ji,jj,jl)
+               ELSE
+                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_low(ji,jj,jl) )
+                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_low(ji,jj,jl) )
+               ENDIF
+            END DO
+         END DO
+
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.  
-               zdiv(ji,jj,jl) =  - (  paa(ji,jj,jl) - paa(ji-1,jj  ,jl)   &
-                  &              +    pbb(ji,jj,jl) - pbb(ji  ,jj-1,jl) )  
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
-
-      DO jl = 1, ipl
-         ! Determine ice masks for before and after tracers 
-         WHERE( pbef(:,:,jl) == 0._wp .AND. paft(:,:,jl) == 0._wp .AND. zdiv(:,:,jl) == 0._wp ) 
-            zmsk(:,:) = 0._wp
-         ELSEWHERE                                                                                   
-            zmsk(:,:) = 1._wp * tmask(:,:,1)
-         END WHERE
-
-         ! 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(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   &
-            &             paft(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  )
-         zbdo(:,:) = MIN( pbef(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   &
-            &             paft(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  )
-
-      z1_dt = 1._wp / pdt
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               !
+               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,jl) ) - MIN( 0., pfu_ho(ji  ,jj,jl) ) &  ! positive/negative part of the flux
+                  & + MAX( 0., pfv_ho(ji,jj-1,jl) ) - MIN( 0., pfv_ho(ji,jj  ,jl) )
+               zneg = MAX( 0., pfu_ho(ji  ,jj,jl) ) - MIN( 0., pfu_ho(ji-1,jj,jl) ) &
+                  & + MAX( 0., pfv_ho(ji,jj  ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) )
+               !
+               IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN
+                  zneg2 = (   pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) &
+                     &    ) * ( 1. - pamsk )
+                  zpos2 = ( - pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) &
+                     &    ) * ( 1. - pamsk )
+               ELSE
+                  zneg2 = 0. ; zpos2 = 0.
+               ENDIF
+               !
+               !                                  ! up & down beta terms
+               IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_low(ji,jj,jl) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt
+               ELSE                         ; zbetup(ji,jj,jl) = 0. ! zbig
+               ENDIF
+               !
+               IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_low(ji,jj,jl) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt
+               ELSE                         ; zbetdo(ji,jj,jl) = 0. ! zbig
+               ENDIF
+               !
+               ! if all the points are outside ice cover
+               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0. ! zbig
+               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 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
             !
-            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  ,jl) ) - MIN( 0., paa(ji  ,jj  ,jl) )   &        ! positive/negative  part of the flux
-                  & + MAX( 0., pbb(ji  ,jj-1,jl) ) - MIN( 0., pbb(ji  ,jj  ,jl) )
-               zneg = MAX( 0., paa(ji  ,jj  ,jl) ) - MIN( 0., paa(ji-1,jj  ,jl) )   &
-                  & + MAX( 0., pbb(ji  ,jj  ,jl) ) - MIN( 0., pbb(ji  ,jj-1,jl) )
-               !
-               zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms
-               zbetup(ji,jj,jl) = ( zup            - paft(ji,jj,jl) ) / ( zpos + zsml ) * zbt
-               zbetdo(ji,jj,jl) = ( paft(ji,jj,jl) - zdo         )    / ( zneg + zsml ) * zbt
             END DO
          END DO
@@ -703,29 +1622,268 @@
       CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
 
-      DO jl = 1, ipl
-         ! monotonic flux in the i & j direction (paa & pbb)
-         ! -------------------------------------
-         DO jj = 2, jpjm1
+      
+      ! monotonic flux in the y direction
+      ! ---------------------------------
+      DO jl = 1, jpl
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
                zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
-               zcu = 0.5  + SIGN( 0.5 , paa(ji,jj,jl) )
+               zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj,jl) )
                !
-               paa(ji,jj,jl) = paa(ji,jj,jl) * ( zcu * zau + ( 1._wp - zcu) * zbu )
-            END DO
-         END DO
-         !
+               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
+
+               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_low(ji,jj,jl)
+
+!!            IF( ji==26 .AND. jj==86) THEN
+!!               WRITE(numout,*) 'coefU',zcoef
+!!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
+!!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * 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,jl) , zbetup(ji,jj+1,jl) )
                zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
-               zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj,jl) )
+               zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj,jl) )
                !
-               pbb(ji,jj,jl) = pbb(ji,jj,jl) * ( zcv * zav + ( 1._wp - zcv) * zbv )
-            END DO
-         END DO
+               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
+
+               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_low(ji,jj,jl)
+
+!!            IF( ji==26 .AND. jj==86) THEN
+!!               WRITE(numout,*) 'coefV',zcoef
+!!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
+!!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * 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,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               ELSE
+                  zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &         ) * tmask(ji,jj,1)
+               ENDIF
+               IF( zzt(ji,jj,jl) < -epsi20 ) THEN
+                  WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj,jl)
+               ENDIF
+            END DO
+         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 (:,:  ), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer
+      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfu_ho       ! high order flux
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux
+      !
+      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
+      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpx       ! tracer slopes 
+      !!----------------------------------------------------------------------
+      !
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond.
+      
+      DO jl = 1, jpl
+         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,jl)
+               Rj  = zslpx(ji  ,jj,jl)
+               Rjp = zslpx(ji+1,jj,jl)
+
+               IF( PRESENT(pfu_ups) ) THEN
+
+                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
+                  ELSE                        ;   Rr = Rjp
+                  ENDIF
+
+                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)     
+                  IF( Rj > 0. ) THEN
+                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  &
+                        &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
+                  ELSE
+                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  &
+                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
+                  ENDIF
+                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + 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,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
+
+               ENDIF
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( 'icedyn_adv_umx', 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 (:,:  ), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer
+      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfv_ho       ! high order flux
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux
+      !
+      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
+      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpy       ! tracer slopes 
+      !!----------------------------------------------------------------------
+      !
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond.
+
+      DO jl = 1, jpl
+         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,jl)
+               Rj  = zslpy(ji,jj  ,jl)
+               Rjp = zslpy(ji,jj+1,jl)
+
+               IF( PRESENT(pfv_ups) ) THEN
+
+                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
+                  ELSE                        ;   Rr = Rjp
+                  ENDIF
+
+                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)     
+                  IF( Rj > 0. ) THEN
+                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  &
+                        &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
+                  ELSE
+                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  &
+                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
+                  ENDIF
+                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + 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,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
+
+               ENDIF
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond.
+      !
+   END SUBROUTINE limiter_y
 
 #else
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rdgrft.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rdgrft.F90	(revision 10419)
@@ -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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg.F90	(revision 10419)
@@ -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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90	(revision 10419)
@@ -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                           
@@ -144,7 +147,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
 !!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zssh_lead_m                     ! array used for the calculation of ice surface slope:
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
       !                                                                 !    ocean surface (ssh_m) if ice is not embedded
-      !                                                                 !    ice top surface if ice is embedded   
+      !                                                                 !    ice bottom surface if ice is embedded   
       REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
@@ -256,13 +259,17 @@
          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
       !
       !------------------------------------------------------------------------------!
       ! 2) Wind / ocean stress, mass terms, coriolis terms
       !------------------------------------------------------------------------------!
-
-      !== embedded sea ice: compute representative ice top surface      ==!
-      !== non-embedded sea ice: use ocean surface for slope calculation ==!
-      zssh_lead_m(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
+      ! sea surface height
+      !    embedded sea ice: compute representative ice top surface
+      !    non-embedded sea ice: use ocean surface for slope calculation
+      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
 
       DO jj = 2, jpjm1
@@ -302,6 +309,6 @@
 
             ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
-            zspgU(ji,jj)    = - zmassU * grav * ( zssh_lead_m(ji+1,jj) - zssh_lead_m(ji,jj) ) * r1_e1u(ji,jj)
-            zspgV(ji,jj)    = - zmassV * grav * ( zssh_lead_m(ji,jj+1) - zssh_lead_m(ji,jj) ) * r1_e2v(ji,jj)
+            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
 
             ! masks
@@ -317,4 +324,43 @@
       CALL lbc_lnk_multi( 'icedyn_rhg_evp', 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( 'icedyn_rhg_evp', 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
@@ -382,7 +428,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
@@ -403,6 +449,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
@@ -452,6 +498,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)
@@ -464,5 +510,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)
@@ -500,6 +546,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)
@@ -512,5 +558,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)
@@ -550,6 +596,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)
@@ -562,5 +608,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)
@@ -598,6 +644,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)
@@ -610,5 +656,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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/iceistate.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/iceistate.F90	(revision 10419)
@@ -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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icestp.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icestp.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icestp.F90	(revision 10419)
@@ -194,13 +194,6 @@
                                         CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
          !
-!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
-!!       moreover it should only be called at the update frequency (as in agrif_ice_update.F90)
-!# if defined key_agrif
-!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid()
-!# endif
                                         CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
-!# if defined key_agrif
-!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()
-!# endif
+         !
          IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs 
          !
@@ -368,9 +361,15 @@
       e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy
       WHERE( a_i_b(:,:,:) >= epsi20 )
-         h_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness
-         h_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness
+         h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:)   ! ice thickness
+         h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:)   ! snw thickness
       ELSEWHERE
          h_i_b(:,:,:) = 0._wp
          h_s_b(:,:,:) = 0._wp
+      END WHERE
+      
+      WHERE( a_ip(:,:,:) >= epsi20 )
+         h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)   ! ice pond thickness
+      ELSEWHERE
+         h_ip_b(:,:,:) = 0._wp
       END WHERE
       !
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icevar.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icevar.F90	(revision 10419)
@@ -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
@@ -953,5 +954,5 @@
    FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
       !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE rhg_evp_rst  ***
+      !!                   ***  ROUTINE ice_var_sshdyn  ***
       !!                     
       !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
@@ -962,5 +963,4 @@
       !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*, 
       !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
-      !!
       !!----------------------------------------------------------------------
       !
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icewri.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icewri.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icewri.F90	(revision 10419)
@@ -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( 'icewri', 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/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/cpl_oasis3.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/cpl_oasis3.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/cpl_oasis3.F90	(revision 10419)
@@ -142,5 +142,5 @@
       INTEGER :: id_part
       INTEGER :: paral(5)       ! OASIS3 box partition
-      INTEGER :: ishape(2,2)    ! shape of arrays passed to PSMILe
+      INTEGER :: ishape(4)    ! shape of arrays passed to PSMILe
       INTEGER :: ji,jc,jm       ! local loop indicees
       CHARACTER(LEN=64) :: zclname
@@ -181,6 +181,8 @@
       !     nl* is set to the global values 1 and jp*glo.
       !
-      ishape(:,1) = (/ 1, nlei-nldi+1 /)
-      ishape(:,2) = (/ 1, nlej-nldj+1 /)
+      ishape(1) = 1
+      ishape(2) = nlei-nldi+1
+      ishape(3) = 1
+      ishape(4) = nlej-nldj+1
       !
       ! ... Allocate memory for data exchange
@@ -242,5 +244,5 @@
 #endif
                   IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
-                  CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 0 /),   &
+                  CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   &
                      &                OASIS_Out          , ishape , OASIS_REAL, nerror )
                   IF ( nerror /= OASIS_Ok ) THEN
@@ -287,5 +289,5 @@
 #endif
                   IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In
-                  CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 0 /),   &
+                  CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   &
                      &                OASIS_In           , ishape , OASIS_REAL, nerror )
                   IF ( nerror /= OASIS_Ok ) THEN
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/sbccpl.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/sbccpl.F90	(revision 10419)
@@ -1667,4 +1667,7 @@
          zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
+      CASE( 'none'      )       ! Not available as for now: needs additional coding below when computing zevap_oce 
+      !                         ! since fields received are not defined with none option
+         CALL ctl_stop( 'STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
       END SELECT
 
@@ -1943,4 +1946,7 @@
             &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       &
             &                     + palbi      (:,:,1) * picefr(:,:) ) )
+      CASE( 'none'      )       ! Not available as for now: needs additional coding  
+      !                         ! since fields received, here zqsr_tot,  are not defined with none option
+         CALL ctl_stop( 'STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_qsr value in namelist namsbc_cpl' )
       END SELECT
       IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/stpctl.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/stpctl.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/stpctl.F90	(revision 10419)
@@ -143,5 +143,5 @@
       !                                   !==  error handling  ==!
       IF( ( ln_ctl .OR. lsomeoce ) .AND. (   &             ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points
-         &  zmax(1) >   15._wp .OR.   &                    ! too large sea surface height ( > 15 m )
+         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m )
          &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s)
          &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity
@@ -161,5 +161,5 @@
          ENDIF
          
-         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 10 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
+         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
          WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2)
          WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3)
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zfechem.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 10419)
@@ -25,5 +25,4 @@
    PUBLIC   p4z_fechem_init   ! called in trcsms_pisces.F90
 
-   LOGICAL          ::   ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
    LOGICAL          ::   ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
    REAL(wp), PUBLIC ::   xlam1        !: scavenging rate of Iron 
@@ -32,6 +31,4 @@
    REAL(wp), PUBLIC ::   kfep         !: rate constant for nanoparticle formation
 
-   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth      !!gm  <<<== DOCTOR names SVP !!!
-
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -47,10 +44,6 @@
       !! ** Purpose :   Compute remineralization/scavenging of iron
       !!
-      !! ** Method  :   2 different chemistry models are available for iron
-      !!                (1) The simple chemistry model of Aumont and Bopp (2006)
-      !!                    based on one ligand and one inorganic form
-      !!                (2) The complex chemistry model of Tagliabue and 
-      !!                    Voelker (2009) based on 2 ligands, 2 inorganic forms
-      !!                    and one particulate form (ln_fechem)
+      !! ** Method  :   A simple chemistry model of iron from Aumont and Bopp (2006)
+      !!                based on one ligand and one inorganic form
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
@@ -74,6 +67,4 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, zFeL1
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zFeL2, zTL2, zFe2, zFeP
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::   zstrn, zstrn2
       !!---------------------------------------------------------------------
       !
@@ -83,12 +74,4 @@
       zFeL1(:,:,:) = 0.
       zTL1 (:,:,:) = 0.
-      IF( ln_fechem ) THEN
-         ALLOCATE( zstrn(jpi,jpj), zstrn2(jpi,jpj) )
-         ALLOCATE( zFe2(jpi,jpj,jpk), zFeL2(jpi,jpj,jpk), zTL2(jpi,jpj,jpk), zFeP(jpi,jpj,jpk) )
-         zFe2 (:,:,:) = 0.
-         zFeL2(:,:,:) = 0.
-         zTL2 (:,:,:) = 0.
-         zFeP (:,:,:) = 0.
-      ENDIF
 
       ! Total ligand concentration : Ligands can be chosen to be constant or variable
@@ -104,134 +87,26 @@
       ENDIF
 
-      IF( ln_fechem ) THEN
-         ! compute the day length depending on latitude and the day
-         zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-         zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-         ! day length in hours
-         zstrn(:,:) = 0.
+      ! ------------------------------------------------------------
+      !  from Aumont and Bopp (2006)
+      ! This model is based on one ligand and Fe' 
+      ! Chemistry is supposed to be fast enough to be at equilibrium
+      ! ------------------------------------------------------------
+      DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-               zargu = MAX( -1., MIN(  1., zargu ) )
-               zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-            END DO
+               zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
+               zkeq            = fekeq(ji,jj,jk)
+               zfesatur        = zTL1(ji,jj,jk) * 1E-9
+               ztfe            = trb(ji,jj,jk,jpfer) 
+               ! Fe' is the root of a 2nd order polynom
+               zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
+                  &              + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
+                  &              + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
+               zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
+               zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
+           END DO
          END DO
-
-         ! Maximum light intensity
-         zstrn2(:,:) = zstrn(:,:) / 24.
-         WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
-         zstrn(:,:) = 24. / zstrn(:,:)
-
-         ! ------------------------------------------------------------
-         ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009)
-         ! This model is based on two ligands, Fe2+, Fe3+ and Fep
-         ! Chemistry is supposed to be fast enough to be at equilibrium
-         ! ------------------------------------------------------------
-         DO jn = 1, 2
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zlight  = etot(ji,jj,jk) * zstrn(ji,jj) * REAL( 2-jn, wp )
-                     zzstrn2 = zstrn2(ji,jj) * REAL( 2-jn, wp ) + (1. - zstrn2(ji,jj) ) * REAL( jn-1, wp )
-                     ! Calculate ligand concentrations : assume 2/3rd of excess goes to
-                     ! strong ligands (L1) and 1/3rd to weak ligands (L2)
-                     ztligand       = ztotlig(ji,jj,jk) - ligand * 1E9
-                     zTL1(ji,jj,jk) =                0.000001 + 0.67 * ztligand
-                     zTL2(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand
-                     ! ionic strength from Millero et al. 1987
-                     zph    = -LOG10( MAX( hi(ji,jj,jk), rtrn) )
-                     zoxy   = trb(ji,jj,jk,jpoxy)
-                     ! Fe2+ oxydation rate from Santana-Casiano et al. (2005)
-                     zkox   = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tempis(ji,jj,jk) + 273.15 )  &
-                     &        - 0.04406 * SQRT( salinprac(ji,jj,jk) ) - 0.002847 * salinprac(ji,jj,jk)
-                     zkox   = ( 10.** zkox ) * spd
-                     zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
-                     ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006)
-                     zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) * (1. - fr_i(ji,jj))
-                     zkph1 = zkph2 / 5.
-                     ! pass the dfe concentration from PISCES
-                     ztfe  = trb(ji,jj,jk,jpfer) * 1e9
-                     ! ----------------------------------------------------------
-                     ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION
-                     ! As shown in Tagliabue and Voelker (2009), Fe3+ is the root of a 3rd order polynom. 
-                     ! ----------------------------------------------------------
-                     ! calculate some parameters
-                     za = 1.0 + ks / kpr
-                     zb = 1.0 + zkph2 / ( zkox )
-                     zc = 1.0 + ( zkph1 + kth ) / ( zkox )
-                     zkappa1 = ( kb1 + zkph1 + kth ) / kl1
-                     zkappa2 = ( kb2 + zkph2 ) / kl2
-                     za2 = zTL2(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za
-                     za1 = zkappa1 * zTL2(ji,jj,jk) * zb / za + zkappa2 * zTL1(ji,jj,jk) * zc / za &
-                      & + zkappa1 * zkappa2 - ( zkappa1 + zkappa2 ) * ztfe / za
-                     za0 = -zkappa1 * zkappa2 * ztfe / za
-                     zp  = za1 - za2 * za2 / 3.
-                     zq  = za2 * za2 * za2 * 2. / 27. - za2 * za1 / 3. + za0
-                     zp3 = zp / 3.
-                     zq2 = zq / 2.
-                     zd  = zp3 * zp3 * zp3 + zq2 * zq2
-                     zr  = zq / ABS( zq ) * SQRT( ABS( zp ) / 3. )
-                     ! compute the roots
-                     IF( zp > 0.) THEN
-                        ! zphi = ASINH( zq / ( 2. * zr * zr * zr ) )
-                        zphi =  zq / ( 2. * zr * zr * zr ) 
-                        zphi = LOG( zphi + SQRT( zphi * zphi + 1 ) )  ! asinh(x) = log(x + sqrt(x^2+1))
-                        zxs  = -2. * zr * SINH( zphi / 3. ) - za1 / 3.
-                     ELSE
-                        IF( zd > 0. ) THEN
-                           zfff = MAX( 1., zq / ( 2. * zr * zr * zr ) )
-                           ! zphi = ACOSH( zfff )
-                           zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) )  ! acosh(x) = log(x + sqrt(x^2-1))
-                           zxs = -2. * zr * COSH( zphi / 3. ) - za1 / 3.
-                        ELSE
-                           zfff = MIN( 1., zq / ( 2. * zr * zr * zr ) )
-                           zphi = ACOS( zfff )
-                           DO jic = 1, 3
-                              zfunc = -2 * zr * COS( zphi / 3. + 2. * REAL( jic - 1, wp ) * rpi / 3. ) - za2 / 3.
-                              IF( zfunc > 0. .AND. zfunc <= ztfe)  zxs = zfunc
-                           END DO
-                        ENDIF
-                     ENDIF
-                     ! solve for the other Fe species
-                     zzFe3  = MAX( 0., zxs )
-                     zzFep  = MAX( 0., ( ks * zzFe3 / kpr ) )
-                     zzFeL1 = MAX( 0., ( zzFe3 * zTL1(ji,jj,jk) ) / ( zkappa1 + zzFe3 ) )
-                     zzFeL2 = (ztfe - za * zzFe3 - zc * zzFeL1 ) / zb
-                     zzFe2  = MAX( 0., ( ( ( zkph1 + kth ) * zzFeL1 + zkph2 * zzFeL2 ) / zkox ) )
-                     zzFep  = ztfe - zzFe3 - zzFe2 - zzFeL1 - zzFeL2
-                     zFe3(ji,jj,jk)  = zFe3(ji,jj,jk)  + zzFe3 * zzstrn2
-                     zFe2(ji,jj,jk)  = zFe2(ji,jj,jk)  + zzFe2 * zzstrn2
-                     zFeL2(ji,jj,jk) = zFeL2(ji,jj,jk) + zzFeL2 * zzstrn2
-                     zFeL1(ji,jj,jk) = zFeL1(ji,jj,jk) + zzFeL1 * zzstrn2
-                     zFeP(ji,jj,jk)  = zFeP(ji,jj,jk)  + zzFeP * zzstrn2
-                  END DO
-               END DO
-            END DO
-         END DO
-      ELSE
-         ! ------------------------------------------------------------
-         ! OLD FE CHEMISTRY ROUTINE from Aumont and Bopp (2006)
-         ! This model is based on one ligand and Fe' 
-         ! Chemistry is supposed to be fast enough to be at equilibrium
-         ! ------------------------------------------------------------
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
-                  zkeq            = fekeq(ji,jj,jk)
-                  zfesatur        = zTL1(ji,jj,jk) * 1E-9
-                  ztfe            = trb(ji,jj,jk,jpfer) 
-                  ! Fe' is the root of a 2nd order polynom
-                  zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
-                     &              + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
-                     &              + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
-                  zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
-                  zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
-              END DO
-            END DO
-         END DO
-         !
-      ENDIF
+      END DO
+         !
 
       zdust = 0.         ! if no dust available
@@ -247,18 +122,13 @@
                &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
                &         + fesol(ji,jj,jk,5) / zhplus )
-               IF( ln_fechem ) THEN
-                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) + zFeP(ji,jj,jk) ) * 1E-9
-                  zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
-                  precip(ji,jj,jk) = 0.0
-               ELSE
-                  zfeequi = zFe3(ji,jj,jk) * 1E-9
-                  zhplus  = max( rtrn, hi(ji,jj,jk) )
-                  fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
+               !
+               zfeequi = zFe3(ji,jj,jk) * 1E-9
+               zhplus  = max( rtrn, hi(ji,jj,jk) )
+               fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
                   &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
                   &         + fesol(ji,jj,jk,5) / zhplus )
-                  zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
-                  ! precipitation of Fe3+, creation of nanoparticles
-                  precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * 1E-9 - fe3sol ) ) * kfep * xstep
-               ENDIF
+               zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
+               ! precipitation of Fe3+, creation of nanoparticles
+               precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * 1E-9 - fe3sol ) ) * kfep * xstep
                !
                ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
@@ -311,7 +181,5 @@
       !  Define the bioavailable fraction of iron
       !  ----------------------------------------
-      IF( ln_fechem ) THEN  ;  biron(:,:,:) = MAX( 0., trb(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
-      ELSE                  ;  biron(:,:,:) = trb(:,:,:,jpfer) 
-      ENDIF
+      biron(:,:,:) = trb(:,:,:,jpfer) 
       !
       IF( ln_ligand ) THEN
@@ -327,5 +195,4 @@
                   zaggliga = zlam1a * xstep * zligco
                   zaggligb = zlam1b * xstep * zligco
-                  tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk)
                   tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb
                   zlcoll3d(ji,jj,jk)  = zaggliga + zaggligb
@@ -334,7 +201,5 @@
          END DO
          !
-         IF( .NOT.ln_fechem) THEN
-            plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trb(:,:,:,jpfer) +rtrn ) ) )
-         ENDIF
+         plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trb(:,:,:,jpfer) +rtrn ) ) )
          !
       ENDIF
@@ -343,19 +208,13 @@
       IF( lk_iomput ) THEN
          IF( knt == nrdttrc ) THEN
-         zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
-         IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
-         IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
-         IF( iom_use("TL1")    )  CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
-         IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
-         IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:)  * 1e9 * tmask(:,:,:) )   ! biron
-         IF( iom_use("FESCAV") )  CALL iom_put("FESCAV" , zscav3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
-         IF( iom_use("FECOLL") )  CALL iom_put("FECOLL" , zcoll3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
-         IF( iom_use("LGWCOLL"))  CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
-         IF( ln_fechem ) THEN
-            IF( iom_use("Fe2")  ) CALL iom_put("Fe2"    , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
-            IF( iom_use("FeL2") ) CALL iom_put("FeL2"   , zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
-            IF( iom_use("FeP")  ) CALL iom_put("FeP"    , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
-            IF( iom_use("TL2")  ) CALL iom_put("TL2"    , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
-         ENDIF
+            zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
+            IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
+            IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
+            IF( iom_use("TL1")    )  CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
+            IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
+            IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:)  * 1e9 * tmask(:,:,:) )   ! biron
+            IF( iom_use("FESCAV") )  CALL iom_put("FESCAV" , zscav3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
+            IF( iom_use("FECOLL") )  CALL iom_put("FECOLL" , zcoll3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
+            IF( iom_use("LGWCOLL"))  CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
          ENDIF
       ENDIF
@@ -367,9 +226,4 @@
       ENDIF
       !
-      IF( ln_fechem )  THEN
-         DEALLOCATE( zstrn, zstrn2 )
-         DEALLOCATE( zFe2, zFeL2, zTL2, zFeP )
-      ENDIF
-      !
       IF( ln_timing )   CALL timing_stop('p4z_fechem')
       !
@@ -391,5 +245,5 @@
       INTEGER ::   ios   ! Local integer 
       !!
-      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand, kfep 
+      NAMELIST/nampisfer/ ln_ligvar, xlam1, xlamdust, ligand, kfep 
       !!----------------------------------------------------------------------
       !
@@ -410,5 +264,4 @@
       IF(lwp) THEN                     ! control print
          WRITE(numout,*) '   Namelist : nampisfer'
-         WRITE(numout,*) '      enable complex iron chemistry scheme      ln_fechem    =', ln_fechem
          WRITE(numout,*) '      variable concentration of ligand          ln_ligvar    =', ln_ligvar
          WRITE(numout,*) '      scavenging rate of Iron                   xlam1        =', xlam1
@@ -418,30 +271,4 @@
       ENDIF
       ! 
-      IF (ln_ligand .AND. ln_fechem) CALL ctl_stop( 'STOP', 'p4z_fechem_init: ln_ligand and ln_fechem are incompatible')
-      !
-      IF( ln_fechem ) THEN             ! set some constants used by the complexe chemistry scheme
-         !
-         spd = 3600. * 24.
-         con = 1.E9
-         ! LIGAND KINETICS (values from Witter et al. 2000)
-         ! Weak (L2) ligands
-         ! Phaeophytin
-         kl2 = 12.2E5  * spd / con
-         kb2 = 12.3E-6 * spd
-         ! Strong (L1) ligands
-         ! Saccharides
-         ! kl1 = 12.2E5  * spd / con
-         ! kb1 = 12.3E-6 * spd
-         ! DFOB-
-         kl1 = 19.6e5  * spd / con
-         kb1 = 1.5e-6  * spd
-         ! pcp and remin of Fe3p
-         ks  = 0.075
-         kpr = 0.05
-         ! thermal reduction of Fe3
-         kth = 0.0048 * 24.
-         !
-      ENDIF
-      !
    END SUBROUTINE p4z_fechem_init
    
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zligand.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 10419)
@@ -25,5 +25,4 @@
    REAL(wp), PUBLIC ::  rlig     !: Remin ligand production
    REAL(wp), PUBLIC ::  prlgw    !: Photochemical of weak ligand
-   REAL(wp), PUBLIC ::  rfep     !: Dissolution rate of FeP
 
    !!----------------------------------------------------------------------
@@ -43,5 +42,5 @@
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zlgwp, zlgwpr, zlgwr, zlablgw, zrfepa, zfepr
+      REAL(wp) ::   zlgwp, zlgwpr, zlgwr, zlablgw
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zligrem, zligpr, zrligprod
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d
@@ -71,16 +70,4 @@
                zligpr(ji,jj,jk)    = zlgwpr
                zrligprod(ji,jj,jk) = zlgwp
-               !
-               ! ----------------------------------------------------------
-               ! Dissolution of nanoparticle Fe
-               ! ----------------------------------------------------------
-               ! dissolution rate is maximal in the presence of light and 
-               ! lower in the aphotici zone
-               ! ! 25 Wm-2 constant
-               zrfepa = rfep * ( 1. - EXP( -1. * etot(ji,jj,jk) / 25. ) ) * (1.- fr_i(ji,jj))
-               zrfepa = MAX( (zrfepa / 10.0), zrfepa ) ! min of 10 days lifetime
-               zfepr  = rfep * xstep * trb(ji,jj,jk,jpfep)
-               tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) - zfepr
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zfepr
                !
             END DO
@@ -130,5 +117,5 @@
       INTEGER ::   ios   ! Local integer 
       !
-      NAMELIST/nampislig/ rlgw, prlgw, rlgs, rfep, rlig
+      NAMELIST/nampislig/ rlgw, prlgw, rlgs, rlig
       !!----------------------------------------------------------------------
       !
@@ -148,5 +135,4 @@
       IF(lwp) THEN                         ! control print
          WRITE(numout,*) '   Namelist : nampislig'
-         WRITE(numout,*) '      Dissolution rate of FeP                      rfep  =', rfep
          WRITE(numout,*) '      Lifetime (years) of weak ligands             rlgw  =', rlgw
          WRITE(numout,*) '      Remin ligand production per unit C           rlig  =', rlig
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zprod.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 10419)
@@ -26,5 +26,4 @@
    PUBLIC   p4z_prod_alloc
 
-   LOGICAL , PUBLIC ::   ln_newprod   !:
    REAL(wp), PUBLIC ::   pislopen     !:
    REAL(wp), PUBLIC ::   pisloped     !:
@@ -156,49 +155,26 @@
       END DO
 
-      IF( ln_newprod ) THEN
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                      ! Computation of production function for Carbon
-                      !  ---------------------------------------------
-                      zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
-                      &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
-                      zpisloped = zpislopeadd(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
-                      &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
-                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
-                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
-                      !  Computation of production function for Chlorophyll
-                      !--------------------------------------------------
-                      zpislopen = zpislopeadn(ji,jj,jk) / ( zprmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-                      zpisloped = zpislopeadd(ji,jj,jk) / ( zprmaxd(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-                      zprnch(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
-                      zprdch(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ELSE
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                      ! Computation of production function for Carbon
-                      !  ---------------------------------------------
-                      zpislopen = zpislopeadn(ji,jj,jk)  / ( zprbio(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
-                      zpisloped = zpislopeadd(ji,jj,jk) / ( zprdia(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
-                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
-                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
-                      !  Computation of production function for Chlorophyll
-                      !--------------------------------------------------
-                      zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-                      zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-                      zprnch(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
-                      zprdch(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                   ! Computation of production function for Carbon
+                   !  ---------------------------------------------
+                   zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+                   &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
+                   zpisloped = zpislopeadd(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+                   &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
+                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
+                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
+                   !  Computation of production function for Chlorophyll
+                   !--------------------------------------------------
+                   zpislopen = zpislopeadn(ji,jj,jk) / ( zprmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+                   zpisloped = zpislopeadd(ji,jj,jk) / ( zprmaxd(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+                   zprnch(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
+                   zprdch(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       !  Computation of a proxy of the N/C ratio
@@ -505,5 +481,5 @@
       INTEGER ::   ios   ! Local integer
       !
-      NAMELIST/namp4zprod/ pislopen, pisloped, xadap, ln_newprod, bresp, excretn, excretd,  &
+      NAMELIST/namp4zprod/ pislopen, pisloped, xadap, bresp, excretn, excretd,  &
          &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip
       !!----------------------------------------------------------------------
@@ -525,5 +501,4 @@
       IF(lwp) THEN                         ! control print
          WRITE(numout,*) '   Namelist : namp4zprod'
-         WRITE(numout,*) '      Enable new parame. of production (T/F)   ln_newprod    =', ln_newprod
          WRITE(numout,*) '      mean Si/C ratio                           grosip       =', grosip
          WRITE(numout,*) '      P-I slope                                 pislopen     =', pislopen
@@ -531,8 +506,6 @@
          WRITE(numout,*) '      excretion ratio of nanophytoplankton      excretn      =', excretn
          WRITE(numout,*) '      excretion ratio of diatoms                excretd      =', excretd
-         IF( ln_newprod )  THEN
-            WRITE(numout,*) '      basal respiration in phytoplankton        bresp        =', bresp
-            WRITE(numout,*) '      Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
-         ENDIF
+         WRITE(numout,*) '      basal respiration in phytoplankton        bresp        =', bresp
+         WRITE(numout,*) '      Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
          WRITE(numout,*) '      P-I slope  for diatoms                    pisloped     =', pisloped
          WRITE(numout,*) '      Minimum Chl/C in nanophytoplankton        chlcnm       =', chlcnm
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsbc.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsbc.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsbc.F90	(revision 10419)
@@ -31,5 +31,4 @@
    REAL(wp), PUBLIC ::   dustsolub    !: Solubility of the dust
    REAL(wp), PUBLIC ::   mfrac        !: Mineral Content of the dust
-   REAL(wp), PUBLIC ::   rdustfep     !: Fraction of dust that is dissolvable
    REAL(wp), PUBLIC ::   icefeinput   !: Iron concentration in sea ice
    REAL(wp), PUBLIC ::   wdust        !: Sinking speed of the dust 
@@ -39,6 +38,4 @@
    REAL(wp)         ::   hratio       !: Fe:3He ratio assumed for vent iron supply
    REAL(wp)         ::   distcoast    !: Distance off the coast for Iron from sediments
-   REAL(wp), PUBLIC ::   fep_rats     !: Fep/Fer ratio from sed  sources
-   REAL(wp), PUBLIC ::   fep_rath     !: Fep/Fer ratio from hydro sources
    REAL(wp), PUBLIC ::   lgw_rath     !: Weak ligand ratio from hydro sources
 
@@ -228,5 +225,5 @@
         &                ln_dust, ln_solub, ln_river, ln_ndepo, ln_ironsed, ln_ironice, ln_hydrofe,    &
         &                sedfeinput, distcoast, dustsolub, icefeinput, wdust, mfrac, nitrfix, diazolight, concfediaz, &
-        &                hratio, fep_rats, fep_rath, rdustfep, lgw_rath
+        &                hratio, lgw_rath
       !!----------------------------------------------------------------------
       !
@@ -265,7 +262,4 @@
          WRITE(numout,*) '      Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
          IF( ln_ligand ) THEN
-            WRITE(numout,*) '      Fep/Fer ratio from sed sources            fep_rats   = ', fep_rats
-            WRITE(numout,*) '      Fep/Fer ratio from sed hydro sources      fep_rath   = ', fep_rath
-            WRITE(numout,*) '      Fraction of dust that is dissolvable      rdustfep   = ', rdustfep
             WRITE(numout,*) '      Weak ligand ratio from sed hydro sources  lgw_rath   = ', lgw_rath
          ENDIF
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 10419)
@@ -60,5 +60,4 @@
       REAL(wp) ::  ztrfer, ztrpo4s, ztrdp, zwdust, zmudia, ztemp
       REAL(wp) ::  xdiano3, xdianh4
-      REAL(wp) ::  zwssfep
       !
       CHARACTER (len=25) :: charout
@@ -68,5 +67,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zsidep, zwsfep, zironice
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zsidep, zironice
       !!---------------------------------------------------------------------
       !
@@ -87,5 +86,4 @@
       ALLOCATE( ztrpo4(jpi,jpj,jpk) )
       IF( ln_p5z )    ALLOCATE( ztrdop(jpi,jpj,jpk) )
-      IF( ln_ligand ) ALLOCATE( zwsfep(jpi,jpj) )
 
       zdenit2d(:,:) = 0.e0
@@ -131,11 +129,4 @@
          ELSE
             zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
-         ENDIF
-         IF ( ln_ligand ) THEN
-            IF( ln_solub ) THEN
-               tra(:,:,1,jpfep) = tra(:,:,1,jpfep) + rdustfep * (1.0 - solub(:,:)) * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85
-            ELSE
-               tra(:,:,1,jpfep) = tra(:,:,1,jpfep) + rdustfep * (1.0 - dustsolub) * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85
-            ENDIF
          ENDIF
          zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1 
@@ -216,5 +207,4 @@
             tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
          IF( ln_ligand ) THEN
-            tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + ( hydrofe(:,:,:) * fep_rath ) * rfact2
             tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * lgw_rath ) * rfact2
          ENDIF
@@ -235,14 +225,4 @@
       END DO
       !
-      IF( ln_ligand ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikt  = mbkt(ji,jj)
-               zdep = e3t_n(ji,jj,ikt) / xstep
-               zwsfep(ji,jj)  = MIN( 0.99 * zdep, wsfep(ji,jj,ikt)  )
-            END DO
-         ENDDO
-      ENDIF
-
       IF( .NOT.lk_sed ) THEN
 !
@@ -251,5 +231,4 @@
          IF( ln_ironsed ) THEN
                             tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
-            IF( ln_ligand ) tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + ( ironsed(:,:,:) * fep_rats ) * rfact2
             !
             IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
@@ -334,15 +313,4 @@
          END DO
       END DO
-      !
-      IF( ln_ligand ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikt     = mbkt(ji,jj)
-               zdep    = xstep / e3t_n(ji,jj,ikt) 
-               zwssfep = zwsfep(ji,jj) * zdep
-               tra(ji,jj,ikt,jpfep) = tra(ji,jj,ikt,jpfep) - trb(ji,jj,ikt,jpfep) * zwssfep
-            END DO
-         END DO
-      ENDIF
       !
       IF( ln_p5z ) THEN
@@ -524,5 +492,4 @@
       !
       IF( ln_p5z )    DEALLOCATE( ztrpo4, ztrdop )
-      IF( ln_ligand ) DEALLOCATE( zwsfep )
       !
       IF( ln_timing )  CALL timing_stop('p4z_sed')
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsink.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 10419)
@@ -35,5 +35,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer            !: Small BFe sinking fluxes
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer2           !: Big iron sinking fluxes
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfep      !: Fep sinking fluxes
 
    INTEGER  :: ik100
@@ -126,11 +125,4 @@
       ENDIF
 
-      IF( ln_ligand ) THEN
-         wsfep (:,:,:) = wfep
-         !
-         sinkfep(:,:,:) = 0.e0
-         CALL trc_sink( kt, wsfep, sinkfep , jpfep, rfact2 )
-      ENDIF
-
      ! Total carbon export per year
      IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
@@ -214,5 +206,5 @@
       !!                     ***  ROUTINE p4z_sink_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(3)
+      INTEGER :: ierr(2)
       !!----------------------------------------------------------------------
       !
@@ -224,8 +216,6 @@
          &      sinkfer(jpi,jpj,jpk)                                            , STAT=ierr(1) )                
          !
-      IF( ln_ligand ) ALLOCATE( sinkfep(jpi,jpj,jpk)                            , STAT=ierr(2) )  
-         
       IF( ln_p5z    ) ALLOCATE( sinkingn(jpi,jpj,jpk), sinking2n(jpi,jpj,jpk)   ,     &
-         &                      sinkingp(jpi,jpj,jpk), sinking2p(jpi,jpj,jpk)   , STAT=ierr(3) )
+         &                      sinkingp(jpi,jpj,jpk), sinking2p(jpi,jpj,jpk)   , STAT=ierr(2) )
       !
       p4z_sink_alloc = MAXVAL( ierr )
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10419)
@@ -74,4 +74,5 @@
             CALL p4z_che                              ! initialize the chemical constants
             CALL ahini_for_at(hi)   !  set PH at kt=nit000
+            t_oce_co2_flx_cum = 0._wp
         ELSE
             CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields
@@ -188,6 +189,5 @@
       !!
       NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
-         &                   wfep, ldocp, ldocz, lthet,  &
-         &                   no3rat3, po4rat3
+         &                   ldocp, ldocz, lthet, no3rat3, po4rat3
          !
       NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
@@ -223,5 +223,4 @@
          WRITE(numout,*) '      Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
          IF( ln_ligand ) THEN
-            WRITE(numout,*) '      FeP sinking speed                              wfep   =', wfep
             IF( ln_p4z ) THEN
                WRITE(numout,*) '      Phyto ligand production per unit doc           ldocp  =', ldocp
@@ -493,5 +492,4 @@
             &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
             &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3    
-         IF( ln_ligand)  zwork(:,:,:) = zwork(:,:,:) + trn(:,:,:,jpfep)                
          !
          ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/par_pisces.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/par_pisces.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/par_pisces.F90	(revision 10419)
@@ -56,5 +56,4 @@
    INTEGER, PUBLIC ::   jpgop     !: Big phosphorus particles Concentration
    INTEGER, PUBLIC ::   jplgw     !: Weak Ligands
-   INTEGER, PUBLIC ::   jpfep     !: Fe nanoparticle
 
    !!---------------------------------------------------------------------
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/sms_pisces.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/sms_pisces.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/sms_pisces.F90	(revision 10419)
@@ -47,5 +47,4 @@
    REAL(wp) ::   xkmort            !: ???
    REAL(wp) ::   ferat3            !: ???
-   REAL(wp) ::   wfep              !: ???
    REAL(wp) ::   ldocp             !: ???
    REAL(wp) ::   ldocz             !: ???
@@ -90,7 +89,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio3   !: POC sinking speed 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio4   !: GOC sinking speed
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsfep
-
-
 
    !!*  SMS for the organic matter
@@ -175,7 +171,6 @@
          ! 
          IF( ln_ligand ) THEN
-           ALLOCATE( plig(jpi,jpj,jpk)  , wsfep(jpi,jpj,jpk)  ,   STAT=ierr(8) )
+           ALLOCATE( plig(jpi,jpj,jpk)  ,                         STAT=ierr(8) )
          ENDIF
-         !
       ENDIF
       !
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/trcini_pisces.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/trcini_pisces.F90	(revision 10417)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/trcini_pisces.F90	(revision 10419)
@@ -168,5 +168,4 @@
         IF( cltra == 'PFe'      )   jppfe = jn      !: Picophytoplankton Fe biomass
         IF( cltra == 'LGW'      )   jplgw = jn      !: Weak ligands
-        IF( cltra == 'LFe'      )   jpfep = jn      !: Fe nanoparticle
       END DO
 
@@ -216,5 +215,4 @@
          IF( ln_ligand) THEN
             trn(:,:,:,jplgw) = 0.6E-9
-            trn(:,:,:,jpfep) = 0. * 5.e-6
          ENDIF
          IF( ln_p5z ) THEN
