Index: NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/namelist_ice_ref
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/namelist_ice_ref	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/namelist_ice_ref	(revision 10399)
@@ -49,7 +49,8 @@
 &namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
-   ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
-   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection)
-   ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection)
+   ln_dynALL        = .true.          !  dyn.: full ice dynamics                  (rheology + advection + ridging/rafting + correction)
+   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections     (rheology + advection)
+   ln_dynADV1D      = .false.         !  dyn.: only advection 1D                  (Schar & Smolarkiewicz 1996 test case)
+   ln_dynADV2D      = .false.         !  dyn.: only advection 2D w prescribed vel.(rn_uvice + advection)
       rn_uice       =   0.5           !        prescribed ice u-velocity
       rn_vice       =   0.5           !        prescribed ice v-velocity
Index: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/ice.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/ice.F90	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/ice.F90	(revision 10399)
@@ -145,4 +145,8 @@
    INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
+   !
+   !                                     !!** ice-advection namelist (namdyn_adv) **
+   LOGICAL , PUBLIC ::   ln_adv_Pra       !: Prather        advection scheme
+   LOGICAL , PUBLIC ::   ln_adv_UMx       !: Ultimate-Macho advection scheme
    !
    !                                     !!** ice-surface forcing namelist (namforcing) **
Index: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90	(revision 10399)
@@ -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
@@ -73,9 +76,11 @@
       INTEGER  ::   ji, jj, jl        ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv
-      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhmax
+      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
       !!--------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('icedyn')
+      ! controls
+      IF( ln_timing    )   CALL timing_start('icedyn')                                                             ! timing
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -88,53 +93,54 @@
          tau_icebfr(:,:) = 0._wp
          DO jl = 1, jpl
-            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
-         END DO
-      ENDIF
-      !
-      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 = 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) ) )
+            DO ji = fs_2, fs_jpim1
+               zhi_max(ji,jj,jl) = MAX( epsi20, h_i_b(ji,jj,jl), h_i_b(ji+1,jj  ,jl), h_i_b(ji  ,jj+1,jl), &
+                  &                                              h_i_b(ji-1,jj  ,jl), h_i_b(ji  ,jj-1,jl), &
+                  &                                              h_i_b(ji+1,jj+1,jl), h_i_b(ji-1,jj-1,jl), &
+                  &                                              h_i_b(ji+1,jj-1,jl), h_i_b(ji-1,jj+1,jl) )
+               zhs_max(ji,jj,jl) = MAX( epsi20, h_s_b(ji,jj,jl), h_s_b(ji+1,jj  ,jl), h_s_b(ji  ,jj+1,jl), &
+                  &                                              h_s_b(ji-1,jj  ,jl), h_s_b(ji  ,jj-1,jl), &
+                  &                                              h_s_b(ji+1,jj+1,jl), h_s_b(ji-1,jj-1,jl), &
+                  &                                              h_s_b(ji+1,jj-1,jl), h_s_b(ji-1,jj+1,jl) )
             END DO
          END DO
       END DO
-      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. )
+      CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), '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_dynALL )           !==  all dynamical processes  ==!
+         CALL ice_dyn_rhg   ( kt )                                       ! -- rheology  
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness
+         CALL ice_dyn_rdgrft( kt )                                       ! -- ridging/rafting 
+         CALL ice_cor       ( kt , 1 )                                   ! -- Corrections
 
       CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
-         CALL ice_dyn_rhg   ( kt )                            ! -- rheology  
-         CALL ice_dyn_adv   ( kt )                            ! -- advection of ice
-         CALL Hpiling                                         ! -- simple pile-up (replaces ridging/rafting)
-
-      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities)
+         CALL ice_dyn_rhg   ( kt )                                       ! -- rheology  
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness
+         CALL Hpiling                                                    ! -- simple pile-up (replaces ridging/rafting)
+
+      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
          ALLOCATE( zdivu_i(jpi,jpj) )
-
-         u_ice(:,:) = rn_uice * umask(:,:,1)
-         v_ice(:,:) = rn_vice * vmask(:,:,1)
-         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
-         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
-         ! --- monotonicity test from Lipscomb et al 2004 --- !
+         ! --- 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
+         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
+         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice + correction on ice thickness
 
          ! diagnostics: divergence at T points 
@@ -150,29 +156,55 @@
          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(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
+         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
+         ! ---
+         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max )   ! -- advection of ice + correction on ice thickness
+
+         ! diagnostics: divergence at T points 
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+            END DO
+         END DO
+         CALL lbc_lnk( zdivu_i, 'T', 1. )
+         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
+
+         DEALLOCATE( zdivu_i )
+
       END SELECT
-      !
-      IF( ln_timing )   CALL timing_stop('icedyn')
+       !
+      ! controls
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_timing    )   CALL timing_stop ('icedyn')                                                             ! timing
       !
    END SUBROUTINE ice_dyn
 
 
-   SUBROUTINE Hbig( phmax )
+   SUBROUTINE Hbig( phi_max, phs_max )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE Hbig  ***
       !!
       !! ** Purpose : Thickness correction in case advection scheme creates
-      !!              abnormally tick ice
-      !!
-      !! ** Method  : 1- check whether ice thickness resulting from advection is
-      !!                 larger than the surrounding 9-points before advection
-      !!                 and reduce it if a) divergence or b) convergence & at_i>0.8
-      !!              2- bound ice thickness with hi_max (99m)
+      !!              abnormally tick ice or snow
+      !!
+      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by adapting ice concentration
+      !!              2- check whether snow thickness is larger than the surrounding 9-points
+      !!                 (before advection) and reduce it by sending the excess in the ocean
+      !!              3- check whether snow load deplets the snow-ice interface below sea level$
+      !!                 and reduce it by sending the excess in the ocean
+      !!              4- correct pond fraction to avoid a_ip > a_i
       !!
       !! ** input   : Max thickness of the surrounding 9-points
       !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max   ! max ice thick from surrounding 9-pts
       !
       INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zh, zdv
+      REAL(wp) ::   zhi, zhs, zvs_excess, zfra
       !!-------------------------------------------------------------------
       !
@@ -182,14 +214,43 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax
+               IF ( v_i(ji,jj,jl) > 0._wp ) THEN
                   !
-                  zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl)
-                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)  
-                  !
-                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
-                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
-                     a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
+                  !                               ! -- check h_i -- !
+                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
+                  zhi = v_i (ji,jj,jl) / a_i(ji,jj,jl)
+!!clem                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)  
+!!clem                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
+!!clem                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
+                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
+                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
                   ENDIF
                   !
+                  !                               ! -- check h_s -- !
+                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
+                  zhs = v_s (ji,jj,jl) / a_i(ji,jj,jl)
+                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
+                     zfra = a_i(ji,jj,jl) * phs_max(ji,jj,jl) / MAX( v_s(ji,jj,jl), epsi20 )
+                     !
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     !
+                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra
+                     v_s(ji,jj,jl)   = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
+                  ENDIF           
+                  !
+                  !                               ! -- check snow load -- !
+                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
+                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
+                  !    this imposed mini can artificially make the snow thickness very high (if concentration decreases drastically)
+                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+                  IF( zvs_excess > 0._wp ) THEN
+                     zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 )
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
+                     !
+                     e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra
+                     v_s(ji,jj,jl)   = v_s(ji,jj,jl) - zvs_excess
+                  ENDIF
+                  
                ENDIF
             END DO
@@ -241,6 +302,6 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
-         &             rn_ishlat ,                                             &
+      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
       !!-------------------------------------------------------------------
@@ -259,7 +320,8 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn:'
-         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL      = ', ln_dynFULL
+         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 only         (rn_uvice + adv)                 ln_dynADV       = ', ln_dynADV
+         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
@@ -275,9 +337,11 @@
       ioptio = 0 
       !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
-      IF( ln_dynFULL   ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF
+      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF
       !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
       IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
-      !      !--- advection only with prescribed ice velocities (from namelist)
-      IF( ln_dynADV    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF
+      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
+      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF
+      !      !--- advection 2D only with prescribed ice velocities (from namelist)
+      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF
       !
       IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
Index: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv.F90	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv.F90	(revision 10399)
@@ -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
@@ -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_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90	(revision 10399)
@@ -35,21 +35,49 @@
    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
 
-   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a, z1_a_ups, zua_ups, zva_ups
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_ai, amaxu, amaxv
    
-   ! rachid trick (only for nonosc)
-   LOGICAL :: ll_rachid = .TRUE.
-
+   LOGICAL ll_dens
+
+   ! advect H all the way (and get V=H*A at the end)
+   LOGICAL :: ll_thickness = .FALSE.
+   
+   ! look for 9 points around in nonosc limiter  
+   LOGICAL :: ll_9points = .FALSE.  ! false=better h?
+
+   ! use HgradU in nonosc limiter  
+   LOGICAL :: ll_HgradU = .TRUE.   ! no effect?
+
+   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
+   LOGICAL :: ll_neg = .TRUE.      ! keep TRUE
+   
+   ! limit the fluxes
+   LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !!
+   LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D
+   LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D
+   LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D
+
+   ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc)
+   LOGICAL :: ll_clem   = .TRUE.   ! simpler than rachid and works
+
+   ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*)
+   LOGICAL :: ll_gurvan = .FALSE.  ! must be false for 1D case !!
+
+   ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-)
+   LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h
+
+   ! advect (or not) open water. If not, retrieve it from advection of A
+   LOGICAL :: ll_ADVopw = .FALSE.  !
+   
    ! alternate directions for upstream
-   ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false"
-   ! clem: needs to be set to true in 2D when using prelimiter (otherwise "wavy solutions" are created)
    LOGICAL :: ll_upsxy = .TRUE.
+
+   ! alternate directions for high order
+   LOGICAL :: ll_hoxy = .TRUE.
    
-   ! prelimiter
-   ! clem: use it to avoid overshoot in H
-   LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => better than Devore especially in 2D
-   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11
+   ! 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)
-   ! clem: useless
    LOGICAL :: ll_limiter_it2 = .FALSE.
    
@@ -94,11 +122,11 @@
       REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
       REAL(wp) ::   zcfl , zdt
-      REAL(wp) ::   zeps = 0.0_wp           ! shift in concentration to avoid division by 0
-      !                                     !    must be >= 0.01 and the best seems to be 0.1
-      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho, za_ups
-      REAL(wp), DIMENSION(jpi,jpj) ::   zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho
+      REAL(wp), DIMENSION(jpi,jpj) ::   zhvar
+      REAL(wp), DIMENSION(jpi,jpj) ::   zai_b, zai_a, z1_ai_b
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
+      !
       !
       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
@@ -132,8 +160,8 @@
       END DO
 
-      IF(.NOT. ALLOCATED(z1_a)    )   ALLOCATE(z1_a    (jpi,jpj))
-      IF(.NOT. ALLOCATED(z1_a_ups))   ALLOCATE(z1_a_ups(jpi,jpj))
-      IF(.NOT. ALLOCATED(zua_ups) )   ALLOCATE(zua_ups (jpi,jpj))
-      IF(.NOT. ALLOCATED(zva_ups) )   ALLOCATE(zva_ups (jpi,jpj))
+      IF(.NOT. ALLOCATED(z1_ai))       ALLOCATE(z1_ai(jpi,jpj))
+      IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj))
+      IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj))
+
       !---------------!
       !== advection ==!
@@ -141,86 +169,101 @@
       DO jt = 1, icycle
 
-         zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:) 
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:), zua_ups, zva_ups, za_ups, zua_ho, zva_ho )   ! Open water area 
-
+         IF( ll_ADVopw ) THEN
+            ll_dens=.FALSE.
+            zamsk = 1._wp
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area 
+         ELSE
+            zai_b(:,:) = SUM( pa_i(:,:,:), dim=3 )
+         ENDIF
+         
          DO jl = 1, jpl
-            ! to avoid a problem with the limiter nonosc when A gets close to 0
-            pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1)
-            !
-            WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl)
-            ELSEWHERE                        ;   z1_a(:,:) = 0.
+            !
+            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_i(:,:,jl)
+            ELSEWHERE                         ;   z1_ai_b(:,:) = 0.
             END WHERE
             !
-            zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:) 
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! Ice area
-
-            ! 1/A_ups
-!!            IF( .NOT. ll_rachid ) THEN
-               WHERE( za_ups(:,:) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / za_ups(:,:)
-               ELSEWHERE                       ;   z1_a_ups(:,:) = 0.
-               END WHERE
-!!            ELSE
-!!               WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / pa_i(:,:,jl)
-!!               ELSEWHERE                        ;   z1_a_ups(:,:) = 0.
-!!               END WHERE              
-!!            ENDIF
-!!
-!!            IF( ll_rachid )   zua_ups = zua_ho ; zva_ups = zva_ho
-            zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl), zua_ups, zva_ups )              ! Ice volume
-            zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl), zua_ups, zva_ups )              ! Snw volume
-            zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl), zua_ups, zva_ups )              ! Salt content
-            zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl), zua_ups, zva_ups )              ! Age content
+            IF( ll_zeroup2 ) THEN
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), &
+                        &                              pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) )
+                     amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), &
+                        &                              pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) )
+                 END DO
+               END DO
+               CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.)
+            ENDIF
+            !
+            zamsk = 1._wp
+            ll_dens=.TRUE.
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ho, zva_ho ) ! Ice area
+            ll_dens=.FALSE.
+
+            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_i(:,:,jl)
+            ELSEWHERE                         ;   z1_ai(:,:) = 0.
+            END WHERE              
+
+            IF( ll_thickness ) THEN
+               zua_ho(:,:) = zudy(:,:)
+               zva_ho(:,:) = zvdx(:,:)
+            ENDIF
+            
+            zamsk = 0._wp ; zhvar(:,:) = pv_i (:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_i (:,:,jl) )              ! Ice volume
+            IF( ll_thickness )   pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
+            
+            zamsk = 0._wp ; zhvar(:,:) = pv_s (:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_s (:,:,jl) )              ! Snw volume
+            IF( ll_thickness )   pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
+
+            zamsk = 0._wp ; zhvar(:,:) = psv_i(:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), psv_i(:,:,jl) )              ! Salt content
+
+            zamsk = 0._wp ; zhvar(:,:) = poa_i(:,:,jl) * z1_ai_b(:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), poa_i(:,:,jl) )              ! Age content
+
             DO jk = 1, nlay_i
-               zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:)
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl), zua_ups, zva_ups )           ! Ice heat content
-            END DO
+               zamsk = 0._wp ; zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai_b(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_i(:,:,jk,jl) )           ! Ice heat content
+            END DO
+
             DO jk = 1, nlay_s
-               zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:)
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl), zua_ups, zva_ups )           ! Snw heat content
-            END DO
-            !
-            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_a...)
-               ! to avoid a problem with the limiter nonosc when A gets close to 0
-               pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1)
+               zamsk = 0._wp ; zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai_b(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content
+            END DO
+            !
+            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_ai_b...)
                !
-               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_ip(:,:,jl)
-               ELSEWHERE                         ;   z1_a(:,:) = 0.
+               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai_b(:,:) = 1._wp / pa_ip(:,:,jl)
+               ELSEWHERE                          ;   z1_ai_b(:,:) = 0.
                END WHERE
                !
-               zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:) 
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! mp fraction
-
-               ! 1/A_ups
-               WHERE( za_ups(:,:) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / za_ups(:,:)
-               ELSEWHERE                       ;   z1_a_ups(:,:) = 0.
-               END WHERE
-
-!!            IF( ll_rachid )   zua_ups = zua_ho ; zva_ups = zva_ho
-               zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:)
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl), zua_ups, zva_ups )              ! mp volume
-            ENDIF
-            !
-            ! to avoid a problem with the limiter nonosc when A gets close to 0
+               zamsk = 1._wp
+               ll_dens=.TRUE.
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ho, zva_ho ) ! mp fractio!n
+               ll_dens=.FALSE.
+
+               WHERE( pa_ip(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_ip(:,:,jl)
+               ELSEWHERE                          ;   z1_ai(:,:) = 0.
+               END WHERE              
+
+               zamsk = 0._wp ; zhvar(:,:) = pv_ip(:,:,jl) * z1_ai_b(:,:)
+               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_ip(:,:,jl) )              ! mp volume
+            ENDIF
+            !
+            !
+         END DO
+
+         IF( .NOT. ll_ADVopw ) THEN
+            zai_a(:,:) = SUM( pa_i(:,:,:), dim=3 )
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1)
-                  pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps &
-                     &             + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
-                     &             ) * tmask(ji,jj,1)
-                  IF ( ln_pnd_H12 ) THEN   ! melt ponds
-                     pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps &
-                        &              + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
-                        &              ) * tmask(ji,jj,1)
-                  ENDIF
-              END DO
-            END DO
-            CALL lbc_lnk_multi( pa_i(:,:,jl), 'T',  1., pa_ip(:,:,jl), 'T',  1. )
-            !
-         END DO
-
+                  pato_i(ji,jj) = pato_i(ji,jj) - ( zai_a(ji,jj) - zai_b(ji,jj) ) &
+                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
+               END DO
+            END DO
+            CALL lbc_lnk( pato_i(:,:), 'T',  1. )
+         ENDIF
+         
       END DO
       !
@@ -228,5 +271,5 @@
 
    
-   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ups, pva_ups, pa_ups, pua_ho, pva_ho )
+   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -249,8 +292,6 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt             ! tracer field
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt             ! tracer field
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pua_ups, pva_ups ! upstream u*a fluxes or u
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pa_ups         ! concentration advected upstream
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
       !
@@ -264,5 +305,17 @@
       !  upstream (_ups) advection with initial mass fluxes
       ! ---------------------------------------------------
-      IF( ll_rachid )    zfu_ups=0.; zfv_ups=0.
+      IF( ll_clem )    zfu_ups=0.; zfv_ups=0.
+
+      IF( ll_gurvan .AND. pamsk==0. ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) &
+                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+            END DO
+         END DO
+         CALL lbc_lnk( pt, 'T', 1. )
+      ENDIF
+
+      
       IF( .NOT. ll_upsxy ) THEN
 
@@ -270,6 +323,11 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1
-               zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj)
-               zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1)
+               IF( ll_clem ) THEN
+                  zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
+               ELSE
+                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
+               ENDIF
             END DO
          END DO
@@ -277,6 +335,6 @@
       ELSE
          ! 1 if advection of A
-         ! z1_a_ups already defined IF advection of other tracers
-         IF( pamsk == 1. )   z1_a_ups(:,:) = 1._wp
+         ! z1_ai already defined IF advection of other tracers
+         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp
          !
          IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
@@ -284,12 +342,31 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  IF( ll_clem ) THEN
+                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  ELSE
+                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
+                  ENDIF
                END DO
             END DO
+            
             ! first guess of tracer content from u-flux
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) &
-                     &         * tmask(ji,jj,1) * z1_a_ups(ji,jj)
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
+                           &         * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) &
+                        &         * tmask(ji,jj,1)
+                  ENDIF
+                  IF( ji==26 .AND. jj==86) THEN
+                     WRITE(numout,*) '************************'
+                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -299,21 +376,45 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj+1)
+                  IF( ll_clem ) THEN
+                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1)
+                  ELSE
+                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
+                  ENDIF
                END DO
             END DO
-            !
+
+         !
          ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
             ! flux in y-direction
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1)
+                  IF( ll_clem ) THEN
+                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
+                  ELSE
+                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
+                  ENDIF
                END DO
             END DO
+
             ! first guess of tracer content from v-flux
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
-                     &         * tmask(ji,jj,1) * z1_a_ups(ji,jj) 
-               END DO
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
+                        &                        + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
+                        &         * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
+                        &         * tmask(ji,jj,1)
+                  ENDIF
+                  IF( ji==26 .AND. jj==86) THEN
+                     WRITE(numout,*) '************************'
+                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
+                  ENDIF
+                END DO
             END DO
             CALL lbc_lnk( zpt, 'T', 1. )
@@ -322,5 +423,9 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * zpt(ji+1,jj)
+                  IF( ll_clem ) THEN
+                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj)
+                  ELSE
+                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
+                  ENDIF
                END DO
             END DO
@@ -330,20 +435,14 @@
       ENDIF
 
-      ! Rachid trick
-      ! ------------
-      IF( ll_rachid .AND. kn_limiter /= 1 )  &
-         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_rachid incompatible with limiters other than nonosc' )
-      IF( ll_rachid ) THEN
-         call lbc_lnk (zfu_ups,'U',-1.)
-         call lbc_lnk (zfv_ups,'V',-1.)
-      IF( pamsk == 0. ) THEN
-         WHERE( ABS( puc(:,:) ) > 0._wp )   ;   zfu_ups(:,:) = zfu_ups(:,:) / puc(:,:) * pu(:,:)
-         ELSEWHERE                          ;   zfu_ups(:,:) = 0._wp
-         END WHERE
-
-         WHERE( ABS( pvc(:,:) ) > 0._wp )   ;   zfv_ups(:,:) = zfv_ups(:,:) / pvc(:,:) * pv(:,:)
-         ELSEWHERE                          ;   zfv_ups(:,:) = 0._wp
-         END WHERE         
-      ENDIF
+      IF( ll_clem .AND. kn_limiter /= 1 )  &
+         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' )
+
+      IF( ll_zeroup2 ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( amaxu(ji,jj) == 0._wp )   zfu_ups(ji,jj) = 0._wp
+               IF( amaxv(ji,jj) == 0._wp )   zfv_ups(ji,jj) = 0._wp
+            END DO
+         END DO
       ENDIF
 
@@ -353,8 +452,17 @@
             ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
                &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj)
-            IF( ll_rachid ) THEN
-               zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+            IF( ll_clem ) THEN
+               IF( ll_gurvan ) THEN
+                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+               ELSE
+                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                                      + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
             ELSE
                zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+            ENDIF
+            IF( ji==26 .AND. jj==86) THEN
+               WRITE(numout,*) '**************************'
+               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj)
             ENDIF
          END DO
@@ -377,15 +485,81 @@
          !
       END SELECT
-         
+
+      IF( ll_thickness ) THEN
+         ! final trend with corrected fluxes
+         ! ------------------------------------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               IF( ll_gurvan ) THEN
+                  ztra       = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 
+               ELSE
+                  ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  & 
+                     &           + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) &
+                     &           + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj)
+               ENDIF
+               pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+
+               IF( pt(ji,jj) < -epsi20 ) THEN
+                  WRITE(numout,*) 'T<0 ',pt(ji,jj)
+               ENDIF
+               
+               IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 )   pt(ji,jj) = 0._wp
+               
+               IF( ji==26 .AND. jj==86) THEN
+                  WRITE(numout,*) 'zt high order',pt(ji,jj)
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk( pt, 'T',  1. )
+      ENDIF
+   
+      ! Rachid trick
+      ! ------------
+     IF( ll_clem ) THEN
+         IF( pamsk == 0. ) THEN
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
+                     zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj)
+                     zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj)
+                  ELSE
+                     zfu_ho (ji,jj) = 0._wp
+                     zfu_ups(ji,jj) = 0._wp
+                  ENDIF
+                  !
+                  IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
+                     zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj)
+                     zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj)
+                  ELSE
+                     zfv_ho (ji,jj) = 0._wp  
+                     zfv_ups(ji,jj) = 0._wp  
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDIF
+      ENDIF
+
+      IF( ll_zeroup5 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, fs_jpim1   ! vector opt.
+               zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                      - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1)
+               IF( zpt(ji,jj) < 0. ) THEN
+                  zfu_ho(ji,jj)   = zfu_ups(ji,jj)
+                  zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj)
+                  zfv_ho(ji,jj)   = zfv_ups(ji,jj)
+                  zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1)
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. )
+      ENDIF
+
       ! output upstream trend of concentration and high order fluxes
       ! ------------------------------------------------------------
-      IF( pamsk == 1. ) THEN
-         ! upstream trend of concentration
-         pa_ups(:,:) = zt_ups(:,:)
-         ! upstream and high order u*a
+      IF( ll_dens ) THEN
+         ! high order u*a
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1
-               pua_ups(ji,jj) = zfu_ups(ji,jj)
-               pva_ups(ji,jj) = zfv_ups(ji,jj)
                pua_ho (ji,jj) = zfu_ho (ji,jj)
                pva_ho (ji,jj) = zfv_ho (ji,jj)
@@ -395,15 +569,34 @@
          !!CALL lbc_lnk( pva_ho, 'V',  -1. )
       ENDIF
+
+
+      IF( .NOT.ll_thickness ) THEN
+         ! final trend with corrected fluxes
+         ! ------------------------------------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1 
+               ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
+                  &     * r1_e1e2t(ji,jj) * pdt  
+
+               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN
+               !!   ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
+               !!      &     * r1_e1e2t(ji,jj) * pdt                 
+               !!ENDIF
+               !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN
+               !!   WRITE(numout,*) 'Tc<0 ',ptc(ji,jj)+ztra
+               !!   ztra = 0._wp
+               !!ENDIF
+
+               ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1)
+                             
+               IF( ji==26 .AND. jj==86) THEN
+                  WRITE(numout,*) 'ztc high order',ptc(ji,jj)
+               ENDIF
+               
+            END DO
+         END DO
+         CALL lbc_lnk( ptc, 'T',  1. )
+      ENDIF
       
-      ! final trend with corrected fluxes
-      ! ------------------------------------
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1 
-            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
-               &         ) * r1_e1e2t(ji,jj)  
-            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      CALL lbc_lnk( ptc, 'T',  1. )
       !
    END SUBROUTINE adv_umx
@@ -448,5 +641,9 @@
          END DO
          IF    ( kn_limiter == 1 ) THEN
-            CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+            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 )
@@ -459,4 +656,6 @@
       ELSE                    !-- alternate directions --!
          !
+         IF( pamsk == 1. )   z1_ai(:,:) = 1._wp
+         !
          IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
             !
@@ -473,5 +672,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
+                  IF( ll_clem ) THEN
+                     zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1)
+                  ELSE                     
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -481,5 +684,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
+                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
                END DO
             END DO
@@ -501,5 +704,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
+                  IF( ll_clem ) THEN
+                     zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1)
+                  ELSE
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -509,5 +716,5 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
+                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
                END DO
             END DO
@@ -516,5 +723,9 @@
 
          ENDIF
-         IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         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
          
       ENDIF
@@ -564,17 +775,51 @@
          IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
          !                                                        !--  advective form update in zzt  --!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
-                  &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
-            END DO
-         END DO
-         CALL lbc_lnk( zzt, 'T', 1. )
+
+         IF( ll_1stguess_clem ) THEN
+
+            ! first guess of tracer content from u-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+
+         ELSE
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_gurvan ) THEN
+                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj)
+                  ELSE
+                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
+                  ENDIF
+                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+         ENDIF
+         !
          !                                                        !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
+         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  ==!
@@ -586,27 +831,71 @@
          IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
          !                                                        !--  advective form update in zzt  --!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
-                  &                    - pt  (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
-            END DO
-         END DO
-         CALL lbc_lnk( zzt, 'T', 1. )
+         IF( ll_1stguess_clem ) THEN
+            
+            ! first guess of tracer content from v-flux 
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( ll_clem ) THEN
+                     IF( ll_gurvan ) THEN
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+                     ELSE
+                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
+                           &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+                     ENDIF
+                  ELSE
+                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
+                        &         * tmask(ji,jj,1)
+                  ENDIF
+                END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+            
+         ELSE
+            
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  IF( ll_gurvan ) THEN
+                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj)
+                  ELSE
+                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
+                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
+                  ENDIF
+                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+         ENDIF
+         !
          !                                                        !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
+         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
-            CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+            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)
-            CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
+            IF( ll_clem ) THEN
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
+            ELSE
+               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
+            ENDIF
             ! guess after content field with high order
             DO jj = 2, jpjm1
@@ -618,5 +907,9 @@
             CALL lbc_lnk( zzt, 'T', 1. )
             ! 2nd iteration of nonosc (limit the flux with the limited high order solution)
-            CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
+            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
@@ -678,5 +971,5 @@
       CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
          !        
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
@@ -687,5 +980,5 @@
       CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -697,5 +990,5 @@
       CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -711,5 +1004,5 @@
       CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -725,5 +1018,5 @@
       CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
          !
-         DO jj = 2, jpjm1
+         DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
@@ -742,7 +1035,22 @@
       END SELECT
       !                                                     !-- High order flux in i-direction  --!
+      IF( ll_neg ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               IF( pt_u(ji,jj) < 0._wp ) THEN
+                  pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
+                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
+               ENDIF
+            END DO
+         END DO
+      ENDIF
+
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector opt.
-            pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
+            IF( ll_clem ) THEN
+               pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj)
+            ELSE
+               pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
+            ENDIF
          END DO
       END DO
@@ -806,5 +1114,5 @@
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
+            DO ji = 1, fs_jpim1
                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
                   &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
@@ -814,5 +1122,5 @@
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
+            DO ji = 1, fs_jpim1
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
@@ -824,5 +1132,5 @@
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
+            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)
@@ -837,5 +1145,5 @@
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
+            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)
@@ -850,5 +1158,5 @@
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
          DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
+            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)
@@ -866,7 +1174,22 @@
       END SELECT
       !                                                     !-- High order flux in j-direction  --!
+      IF( ll_neg ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               IF( pt_v(ji,jj) < 0._wp ) THEN
+                  pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
+                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
+               ENDIF
+            END DO
+         END DO
+      ENDIF
+
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector opt.
-            pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
+            IF( ll_clem ) THEN
+               pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj)
+            ELSE
+               pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
+            ENDIF
          END DO
       END DO
@@ -875,5 +1198,5 @@
      
 
-   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
+   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -883,5 +1206,5 @@
       !!
       !! **  Method  :   ... ???
-      !!       warning : ptc and pt_low 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-
@@ -894,29 +1217,90 @@
       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_low      ! before field & upstream guess of after field
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field
       REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux
       REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt   ! local scalars
-      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low
+      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt
       !!----------------------------------------------------------------------
       zbig = 1.e+40_wp
       zsml = epsi20
 
-      ! Rachid trick (upstream already done in macho)
-      ! ------------
-      IF( ll_rachid ) THEN
-      IF( pamsk == 0. ) THEN
-         WHERE( ABS( puc(:,:) ) > 0._wp )   ;   pfu_ho(:,:) = pfu_ho(:,:) / puc(:,:) * pu(:,:)
-         ELSEWHERE                          ;   pfu_ho(:,:) = 0._wp
-         END WHERE
-
-         WHERE( ABS( pvc(:,:) ) > 0._wp )   ;   pfv_ho(:,:) = pfv_ho(:,:) / pvc(:,:) * pv(:,:)
-         ELSEWHERE                          ;   pfv_ho(:,:) = 0._wp
-         END WHERE         
+      IF( ll_zeroup2 ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( amaxu(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp
+               IF( amaxv(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp
+            END DO
+         END DO
       ENDIF
+      
+      IF( ll_zeroup4 ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               IF( pfu_low(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp
+               IF( pfv_low(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp
+            END DO
+         END DO
       ENDIF
+
+
+      IF( ll_zeroup1 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               IF( ll_gurvan ) THEN
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               ELSE
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
+               IF( zzt(ji,jj) < 0._wp ) THEN
+                  pfu_ho(ji,jj)   = pfu_low(ji,jj)
+                  pfv_ho(ji,jj)   = pfv_low(ji,jj)
+                  WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj)
+               ENDIF
+               IF( ji==26 .AND. jj==86) THEN
+                  WRITE(numout,*) 'zzt high order',zzt(ji,jj)
+                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
+                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
+               ENDIF
+               IF( ll_gurvan ) THEN
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               ELSE
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
+               IF( zzt(ji,jj) < 0._wp ) THEN
+                  pfu_ho(ji-1,jj) = pfu_low(ji-1,jj)
+                  pfv_ho(ji,jj-1) = pfv_low(ji,jj-1)
+                  WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj)
+               ENDIF
+               IF( ll_gurvan ) THEN
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+               ELSE
+                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+               ENDIF
+               IF( zzt(ji,jj) < 0._wp ) THEN
+                  WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj)
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )
+      ENDIF
+
       
       ! antidiffusive flux : high order minus low order
@@ -926,5 +1310,5 @@
             pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj)
             pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj)
-         END DO
+          END DO
       END DO
 
@@ -1014,37 +1398,38 @@
       ! Search local extrema
       ! --------------------
-      ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover
+      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover
       DO jj = 1, jpj
          DO ji = 1, jpi
-!!clem            IF    ( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN
-            IF    ( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN
+            IF    ( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN
                zbup(ji,jj) = -zbig
                zbdo(ji,jj) =  zbig
-!!clem            ELSEIF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) /= 0._wp ) THEN
-            ELSEIF( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) > epsi20 ) THEN
+            ELSEIF( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) > 0._wp ) THEN
                zbup(ji,jj) = pt_low(ji,jj)
                zbdo(ji,jj) = pt_low(ji,jj)
-!!clem            ELSEIF( ptc(ji,jj) /= 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN
-            ELSEIF( ptc(ji,jj) > epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN
-               zbup(ji,jj) = ptc(ji,jj)
-               zbdo(ji,jj) = ptc(ji,jj)
+            ELSEIF( pt(ji,jj) > 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN
+               zbup(ji,jj) = pt(ji,jj)
+               zbdo(ji,jj) = pt(ji,jj)
             ELSE
-               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) )
-               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) )
-            ENDIF
-         END DO
-      END DO
-
+               zbup(ji,jj) = MAX( pt(ji,jj) , pt_low(ji,jj) )
+               zbdo(ji,jj) = MIN( pt(ji,jj) , pt_low(ji,jj) )
+            ENDIF
+         END DO
+      END DO
+
+ 
       z1_dt = 1._wp / pdt
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
             !
-!            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), 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) )
-            !
+            IF( .NOT. ll_9points ) THEN
+            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1) )  ! search max/min in neighbourhood
+            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) )
+            !
+            ELSE
             zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1), &  ! search max/min in neighbourhood
                &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  )
             zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), &
                &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  )
+            ENDIF
             !
             zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux
@@ -1053,18 +1438,54 @@
                & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) )
             !
+            IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN
+               zneg2 = (   pt(ji,jj) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
+               zpos2 = ( - pt(ji,jj) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
+            ELSE
+               zneg2 = 0. ; zpos2 = 0.
+            ENDIF
+            !
             !                                  ! up & down beta terms
 !            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
 !            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
-            IF( zpos >= epsi10 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
-            ELSE                      ; zbetup(ji,jj) = 0.
-            ENDIF
-            !
-            IF( zneg >= epsi10 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
-            ELSE                      ; zbetdo(ji,jj) = 0.
-            ENDIF
-            !
-            IF( zbetdo(ji,jj) < 0._wp ) zbetdo(ji,jj)=0.
-            IF( zbetup(ji,jj) < 0._wp ) zbetup(ji,jj)=0.
-            !
+
+            IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt
+            ELSE                         ; zbetup(ji,jj) = 0. ! zbig
+            ENDIF
+            !
+            IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt
+            ELSE                         ; zbetdo(ji,jj) = 0. ! zbig
+            ENDIF
+            !
+            ! if all the points are outside ice cover
+            IF( zup == -zbig )   zbetup(ji,jj) = 0. ! zbig
+            IF( zdo ==  zbig )   zbetdo(ji,jj) = 0. ! zbig            
+            !
+
+            IF( ji==26 .AND. jj==86) THEN
+               WRITE(numout,*) '-----------------'
+               WRITE(numout,*) 'zpos',zpos,zpos2
+               WRITE(numout,*) 'zneg',zneg,zneg2
+               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj)))
+               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj)))
+               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj)))
+               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1)))
+               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt
+               
+               WRITE(numout,*) 'pt',pt(ji,jj)
+               WRITE(numout,*) 'ptim1',pt(ji-1,jj)
+               WRITE(numout,*) 'ptjm1',pt(ji,jj-1)
+               WRITE(numout,*) 'pt_low',pt_low(ji,jj)
+               WRITE(numout,*) 'zbetup',zbetup(ji,jj)
+               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj)
+               WRITE(numout,*) 'zup',zup
+               WRITE(numout,*) 'zdo',zdo
+            ENDIF
             !
          END DO
@@ -1072,25 +1493,4 @@
       CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
 
-      ! Rachid trick
-      ! ------------
-      IF( ll_rachid ) THEN
-      IF( pamsk == 0. ) THEN
-         WHERE( ABS( pu(:,:) ) > 0._wp )
-            pfu_ho (:,:) = pfu_ho (:,:) * puc(:,:) / pu(:,:)
-            pfu_low(:,:) = pfu_low(:,:) * puc(:,:) / pu(:,:)
-         ELSEWHERE
-            pfu_ho (:,:) = 0._wp
-            pfu_low(:,:) = 0._wp
-         END WHERE
-
-         WHERE( ABS( pv(:,:) ) > 0._wp )
-            pfv_ho (:,:) = pfv_ho (:,:) * pvc(:,:) / pv(:,:)
-            pfv_low(:,:) = pfv_low(:,:) * pvc(:,:) / pv(:,:)
-         ELSEWHERE
-            pfv_ho (:,:) = 0._wp
-            pfv_low(:,:) = 0._wp
-         END WHERE         
-      ENDIF
-      ENDIF
       
       ! monotonic flux in the y direction
@@ -1102,5 +1502,14 @@
             zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
             !
-            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj)
+            zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
+            
+            pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj)
+
+            IF( ji==26 .AND. jj==86) THEN
+               WRITE(numout,*) 'coefU',zcoef
+               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
+            ENDIF
+
          END DO
       END DO
@@ -1112,7 +1521,35 @@
             zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
             !
-            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj)
-         END DO
-      END DO
+            zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
+            
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj)
+
+            IF( ji==26 .AND. jj==86) THEN
+               WRITE(numout,*) 'coefV',zcoef
+               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
+               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
+            ENDIF
+         END DO
+      END DO
+
+      ! clem test
+      DO jj = 2, jpjm1
+         DO ji = 2, fs_jpim1   ! vector opt.
+            IF( ll_gurvan ) THEN
+               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+            ELSE
+               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
+                  &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
+                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
+            ENDIF
+            IF( zzt(ji,jj) < -epsi20 ) THEN
+               WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj)
+            ENDIF
+         END DO
+      END DO
+      
+      !
       !
    END SUBROUTINE nonosc_2d
Index: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rdgrft.F90	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rdgrft.F90	(revision 10399)
@@ -192,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_r9947_SI3_advection/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/iceistate.F90	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/iceistate.F90	(revision 10399)
@@ -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_r9947_SI3_advection/tests/ICEADV/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/EXPREF/namelist_ice_cfg	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/EXPREF/namelist_ice_cfg	(revision 10399)
@@ -34,7 +34,8 @@
 &namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
-   ln_dynFULL       = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   ln_dynALL        = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
    ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection)
-   ln_dynADV        = .true.          !  dyn.: only advection w prescribed vel.(rn_uvice + advection)
+   ln_dynADV1D      = .true.          !  dyn.: only advection 1D                  (Schar & Smolarkiewicz 1996 test case)
+   ln_dynADV2D      = .false.         !  dyn.: only advection 2D w prescribed vel.(rn_uvice + advection)
       rn_uice       =   1.            !        prescribed ice u-velocity
       rn_vice       =   0.            !        prescribed ice v-velocity
Index: NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn.F90	(revision 10397)
+++ 	(revision )
@@ -1,313 +1,0 @@
-MODULE icedyn
-   !!======================================================================
-   !!                     ***  MODULE  icedyn  ***
-   !!   Sea-Ice dynamics : master routine for sea ice dynamics 
-   !!======================================================================
-   !! history :  4.0  ! 2018  (C. Rousset)  original code SI3 [aka Sea Ice cube]
-   !!----------------------------------------------------------------------
-#if defined key_si3
-   !!----------------------------------------------------------------------
-   !!   'key_si3'                                       SI3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   ice_dyn       : dynamics of sea ice
-   !!   ice_dyn_init  : initialization and namelist read
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constants
-   USE dom_oce        ! ocean space and time domain
-   USE ice            ! sea-ice: variables
-   USE icedyn_rhg     ! sea-ice: rheology
-   USE icedyn_adv     ! sea-ice: advection
-   USE icedyn_rdgrft  ! sea-ice: ridging/rafting
-   USE icecor         ! sea-ice: corrections
-   USE icevar         ! sea-ice: operations
-   !
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
-   USE lbclnk         ! lateral boundary conditions (or mpp links)
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_dyn        ! called by icestp.F90
-   PUBLIC   ice_dyn_init   ! called by icestp.F90
-   
-   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_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection) 
-   INTEGER, PARAMETER ::   np_dynADV     = 3   ! only advection 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)
-   
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
-   !! $Id: icedyn.F90 8378 2017-07-26 13:55:59Z clem $
-   !! Software governed by the CeCILL licence     (./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_dyn( kt )
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE ice_dyn  ***
-      !!               
-      !! ** Purpose :   this routine manages sea ice dynamics
-      !!
-      !! ** Action : - calculation of friction in case of landfast ice
-      !!             - call ice_dyn_rhg    = rheology
-      !!             - call ice_dyn_adv    = advection
-      !!             - call ice_dyn_rdgrft = ridging/rafting
-      !!             - call ice_cor        = corrections if fields are out of bounds
-      !!--------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! ice time step
-      !!
-      INTEGER  ::   ji, jj, jl        ! dummy loop indices
-      REAL(wp) ::   zcoefu, zcoefv
-      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhmax
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
-      !!--------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('icedyn')
-      !
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*)'ice_dyn: sea-ice dynamics'
-         WRITE(numout,*)'~~~~~~~'
-      ENDIF
-
-      !                      
-      IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction
-         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)
-      DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-!!gm use of MAXVAL here is very probably less efficient than expending the 9 values
-               zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( h_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. )
-      !
-      !
-      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
-
-      CASE ( np_dynFULL )          !==  all dynamical processes  ==!
-         CALL ice_dyn_rhg   ( kt )                            ! -- rheology  
-         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness
-         CALL ice_dyn_rdgrft( kt )                            ! -- ridging/rafting 
-         CALL ice_cor       ( kt , 1 )                        ! -- Corrections
-
-      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
-         CALL ice_dyn_rhg   ( kt )                            ! -- rheology  
-         CALL ice_dyn_adv   ( kt )                            ! -- advection of ice
-         CALL Hpiling                                         ! -- simple pile-up (replaces ridging/rafting)
-
-      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities)
-
-         ALLOCATE( zdivu_i(jpi,jpj) )
-
-         !u_ice(:,:) = rn_uice * umask(:,:,1)
-         !v_ice(:,:) = rn_vice * vmask(:,:,1)
-         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * umask(:,:,1)
-         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * vmask(:,:,1)
-         ! --- monotonicity test from Lipscomb et al 2004 --- !
-         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
-         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-            !   IF( (ji+nimpp-1) < ((jpiglo+1)/6)-5 .OR. (ji+nimpp-1) > (5*(jpiglo+1)/6)+5 ) THEN
-            !      u_ice(ji,jj) = 0._wp
-            !      v_ice(ji,jj) = 0._wp
-            !   ELSE
-                  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)
-             !  ENDIF
-            END DO
-         END DO
-         ! ---
-         CALL ice_dyn_adv   ( kt )                            ! -- advection of ice
-
-         ! diagnostics: divergence at T points 
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( zdivu_i, 'T', 1. )
-         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
-
-         DEALLOCATE( zdivu_i )
-
-      END SELECT
-      !
-      IF( ln_timing )   CALL timing_stop('icedyn')
-      !
-   END SUBROUTINE ice_dyn
-
-
-   SUBROUTINE Hbig( phmax )
-      !!-------------------------------------------------------------------
-      !!                  ***  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)
-      !!
-      !! ** input   : Max thickness of the surrounding 9-points
-      !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts
-      !
-      INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      REAL(wp) ::   zh, zdv
-      !!-------------------------------------------------------------------
-      !
-      CALL ice_var_zapsmall                       !-- zap small areas
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax
-                  !
-                  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)
-                  ENDIF
-                  !
-               ENDIF
-            END DO
-         END DO
-      END DO 
-      !                                           !-- correct pond fraction to avoid a_ip > a_i
-      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
-      !
-   END SUBROUTINE Hbig
-
-
-   SUBROUTINE Hpiling
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE Hpiling  ***
-      !!
-      !! ** Purpose : Simple conservative piling comparable with 1-cat models
-      !!
-      !! ** Method  : pile-up ice when no ridging/rafting
-      !!
-      !! ** input   : a_i
-      !!-------------------------------------------------------------------
-      INTEGER ::   jl         ! dummy loop indices
-      !!-------------------------------------------------------------------
-      !
-      CALL ice_var_zapsmall                       !-- zap small areas
-      !
-      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
-      DO jl = 1, jpl
-         WHERE( at_i(:,:) > epsi20 )
-            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
-         END WHERE
-      END DO
-      !
-   END SUBROUTINE Hpiling
-
-
-   SUBROUTINE ice_dyn_init
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_dyn_init  ***
-      !!
-      !! ** Purpose : Physical constants and parameters linked to the ice
-      !!      dynamics
-      !!
-      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
-      !!       parameter values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist namdyn
-      !!-------------------------------------------------------------------
-      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
-      !!-------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )         ! Namelist namdyn in reference namelist : Ice dynamics
-      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn in reference namelist', lwp )
-      REWIND( numnam_ice_cfg )         ! Namelist namdyn in configuration namelist : Ice dynamics
-      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn in configuration namelist', lwp )
-      IF(lwm) WRITE( numoni, namdyn )
-      !
-      IF(lwp) THEN                     ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
-         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,*)
-      ENDIF
-      !                             !== set the choice of ice dynamics ==!
-      ioptio = 0 
-      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
-      IF( ln_dynFULL   ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   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
-      !
-      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
-      !
-      !                                      !--- Lateral boundary conditions
-      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
-      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
-      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
-      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
-      !
-      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
-      CALL ice_dyn_rhg_init             ! set ice rheology parameters
-      CALL ice_dyn_adv_init             ! set ice advection parameters
-      !
-   END SUBROUTINE ice_dyn_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option         Empty module           NO SI3 sea-ice model
-   !!----------------------------------------------------------------------
-#endif 
-
-   !!======================================================================
-END MODULE icedyn
Index: NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv.F90	(revision 10397)
+++ 	(revision )
@@ -1,198 +1,0 @@
-MODULE icedyn_adv
-   !!======================================================================
-   !!                       ***  MODULE icedyn_adv   ***
-   !!   sea-ice: advection
-   !!======================================================================
-   !! History :  4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
-   !!----------------------------------------------------------------------
-#if defined key_si3
-   !!----------------------------------------------------------------------
-   !!   'key_si3'                                       SI3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   ice_dyn_adv   : advection of sea ice variables
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constant
-   USE dom_oce        ! ocean domain
-   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call
-   USE ice            ! sea-ice: variables
-   USE icevar         ! sea-ice: operations
-   USE icedyn_adv_pra ! sea-ice: advection scheme (Prather)
-   USE icedyn_adv_umx ! sea-ice: advection scheme (ultimate-macho)
-   USE icectl         ! sea-ice: control prints
-   !
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
-   USE timing         ! Timing
-   USE prtctl         ! Print control
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_dyn_adv        ! called by icestp
-   PUBLIC   ice_dyn_adv_init   ! called by icedyn
-
-   INTEGER ::              nice_adv   ! choice of the type of advection scheme
-   !                                        ! associated indices:
-   INTEGER, PARAMETER ::   np_advPRA = 1   ! Prather scheme
-   INTEGER, PARAMETER ::   np_advUMx = 2   ! Ultimate-Macho scheme
-   !
-   ! ** 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   
-   !
-   !! * Substitution
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
-   !! $Id: icedyn_adv.F90 8373 2017-07-25 17:44:54Z clem $
-   !! Software governed by the CeCILL licence     (./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_dyn_adv( kt ) 
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_dyn_adv ***
-      !!                    
-      !! ** purpose : advection of sea ice
-      !!
-      !! ** method  : One can choose between 
-      !!     a) an Ultimate-Macho scheme (with order defined by nn_UMx) => ln_adv_UMx
-      !!     b) and a second order Prather scheme => ln_adv_Pra
-      !!
-      !! ** action :
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! number of iteration
-      !
-      INTEGER ::   jl   ! dummy loop indice
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmask  ! fraction of time step with v_i < 0
-      !!---------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('icedyn_adv')
-      !
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_dyn_adv: sea-ice advection'
-         WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      
-      ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-                           
-      !----------
-      ! Advection
-      !----------
-      SELECT CASE( nice_adv )
-      !                                !-----------------------!
-      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme !
-         !                             !-----------------------!
-         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice,  &
-            &                  ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
-      !                                !-----------------------!
-      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 )
-      END SELECT
-
-      !----------------------------
-      ! 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.
-      !                   These values are usually very small but in some occasions they can also be non-negligible
-      !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix)
-      !
-      ! record the negative values resulting from UMx
-      zmask(:,:) = 0._wp ! keep the init to 0 here
-      DO jl = 1, jpl
-         WHERE( v_i(:,:,jl) < 0._wp )   zmask(:,:) = 1._wp
-      END DO
-      IF( iom_use('iceneg_pres') )   CALL iom_put("iceneg_pres", zmask                                      )  ! fraction of time step with v_i < 0
-      IF( iom_use('iceneg_volu') )   CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 )                )  ! negative ice volume (only)
-      IF( iom_use('iceneg_hfx' ) )   CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. )  &  ! negative ice heat content (only)
-         &                                                                  , dim=4 ), dim=3 ) )* r1_rdtice )  ! -- eq. heat flux [W/m2]
-      !
-      ! ==> conservation is ensured by calling this routine below,
-      !     however the global ice volume is then changed by advection (but errors are very 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 )
-
-      !------------
-      ! diagnostics
-      !------------
-      diag_trp_ei(:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice
-      diag_trp_es(:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice
-      diag_trp_sv(:,:) = SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice
-      diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice
-      diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice
-      IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoi          )   ! ice mass transport
-      IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhos          )   ! snw mass transport
-      IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s)
-      IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                )   ! advected ice heat content (W/m2)
-      IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                )   ! advected snw heat content (W/m2)
-
-      ! controls
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_icectl    )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ')                           ! prints
-      IF( ln_timing    )   CALL timing_stop ('icedyn_adv')                                                             ! timing
-      !
-   END SUBROUTINE ice_dyn_adv
-
-
-   SUBROUTINE ice_dyn_adv_init
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_dyn_adv_init  ***
-      !!
-      !! ** Purpose :   Physical constants and parameters linked to the ice
-      !!                dynamics
-      !!
-      !! ** Method  :   Read the namdyn_adv namelist and check the ice-dynamic
-      !!                parameter values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist namdyn_adv
-      !!-------------------------------------------------------------------
-      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
-      !!
-      NAMELIST/namdyn_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx
-      !!-------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )         ! Namelist namdyn_adv in reference namelist : Ice dynamics
-      READ  ( numnam_ice_ref, namdyn_adv, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_adv in reference namelist', lwp )
-      REWIND( numnam_ice_cfg )         ! Namelist namdyn_adv in configuration namelist : Ice dynamics
-      READ  ( numnam_ice_cfg, namdyn_adv, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_adv in configuration namelist', lwp )
-      IF(lwm) WRITE( numoni, namdyn_adv )
-      !
-      IF(lwp) THEN                     ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_dyn_adv_init: ice parameters for ice dynamics '
-         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namdyn_adv:'
-         WRITE(numout,*) '      type of advection scheme (Prather)             ln_adv_Pra = ', ln_adv_Pra 
-         WRITE(numout,*) '      type of advection scheme (Ulimate-Macho)       ln_adv_UMx = ', ln_adv_UMx 
-         WRITE(numout,*) '         order of the Ultimate-Macho scheme          nn_UMx     = ', nn_UMx
-      ENDIF
-      !
-      !                             !== set the choice of ice advection ==!
-      ioptio = 0 
-      IF( ln_adv_Pra ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advPRA    ;   ENDIF
-      IF( ln_adv_UMx ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advUMx    ;   ENDIF
-      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_adv_init: choose one and only one ice adv. scheme (ln_adv_Pra or ln_adv_UMx)' )
-      !
-      IF( ln_adv_Pra )   CALL adv_pra_init  !* read or initialize all required files
-      !
-   END SUBROUTINE ice_dyn_adv_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option         Empty Module           NO SI3 sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE icedyn_adv
-
Index: NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv_umx.F90	(revision 10397)
+++ 	(revision )
@@ -1,1183 +1,0 @@
-MODULE icedyn_adv_umx
-   !!==============================================================================
-   !!                       ***  MODULE  icedyn_adv_umx  ***
-   !! sea-ice : advection using the ULTIMATE-MACHO scheme
-   !!==============================================================================
-   !! History :  3.6  !  2014-11  (C. Rousset, G. Madec)  Original code
-   !!            4.0  !  2018     (many people)           SI3 [aka Sea Ice cube]
-   !!----------------------------------------------------------------------
-#if defined key_si3
-   !!----------------------------------------------------------------------
-   !!   'key_si3'                                       SI3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme
-   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
-   !!   macho             : ???
-   !!   nonosc            : compute monotonic tracer fluxes by a non-oscillatory algorithm 
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constant
-   USE dom_oce        ! ocean domain
-   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
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
-   USE lbclnk         ! lateral boundary conditions (or mpp links)
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
-      
-   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
-   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
-
-   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a
-   
-   ! alternate directions for upstream
-   ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false"
-   LOGICAL :: ll_upsxy = .FALSE.
-   
-   ! prelimiter
-   LOGICAL :: ll_prelimiter = .TRUE.
-   LOGICAL :: ll_prelimiter_zalesak = .TRUE.   ! from: Zalesak(1979) eq. 14 => better than Devore
-   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11
-
-   ! iterate on the limiter (only nonosc)
-   LOGICAL :: ll_limiter_it2 = .TRUE.
-   
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
-   !! $Id: icedyn_adv_umx.F90 4499 2014-02-18 15:14:31Z timgraham $
-   !! Software governed by the CeCILL licence     (./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   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  ***
-      !! 
-      !! **  Purpose :   Compute the now trend due to total advection of 
-      !!                 tracers and add it to the general trend of tracer equations
-      !!                 using an "Ultimate-Macho" scheme
-      !!
-      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
-      !!----------------------------------------------------------------------
-      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
-      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
-      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
-      !
-      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
-      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
-      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
-      REAL(wp) ::   zcfl , zdt
-      REAL(wp) ::   zeps = 0.0_wp           ! shift in concentration to avoid division by 0
-      !                                     !    must be >= 0.01 and the best seems to be 0.1
-      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zfu, zfv
-      REAL(wp), DIMENSION(jpi,jpj) ::   zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip 
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
-      !
-      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
-      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
-      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
-      IF( lk_mpp )   CALL mpp_max( zcfl )
-
-      IF( zcfl > 0.5 ) THEN   ;   icycle = 2 
-      ELSE                    ;   icycle = 1 
-      ENDIF
-
-      zdt = rdt_ice / REAL(icycle)
-
-      ! --- transport --- !
-      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
-      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
-
-      ! --- define velocity for advection: u*grad(H) --- !
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
-            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
-            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
-            ENDIF
-
-            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
-            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
-            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
-            ENDIF
-         END DO
-      END DO
-
-      IF(.NOT. ALLOCATED(z1_a )) ALLOCATE(z1_a (jpi,jpj))
-      !---------------!
-      !== advection ==!
-      !---------------!
-      DO jt = 1, icycle
-
-         zamsk = 1._wp
-         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area 
-
-         DO jl = 1, jpl
-            ! to avoid a problem with the limiter nonosc when A gets close to 0
-            pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1)
-            !
-            WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl)
-            ELSEWHERE                        ;   z1_a(:,:) = 0.   !!;   pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0. 
-            END WHERE
-            !
-            zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zfu, zfv ) ! Ice area
-            !!zfu = zudy ; zfv = zvdx
-            !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) )
-            zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl) )              ! Ice volume
-            zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) )              ! Snw volume
-            zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) )              ! Salt content
-            zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:)
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl) )              ! Age content
-            DO jk = 1, nlay_i
-               zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:)
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl) )           ! Ice heat content
-            END DO
-            DO jk = 1, nlay_s
-               zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:)
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl) )           ! Snw heat content
-            END DO
-            !
-            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_a...)
-               ! to avoid a problem with the limiter nonosc when A gets close to 0
-               pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1)
-               !
-               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_ip(:,:,jl)
-               ELSEWHERE                         ;   z1_a(:,:) = 0.
-               END WHERE
-               !
-               zamsk = 1._wp
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zfu, zfv ) ! mp fraction
-               zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:)
-               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl) )              ! mp volume
-            ENDIF
-            !
-            ! to avoid a problem with the limiter nonosc when A gets close to 0
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1)
-                  pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps &
-                     &             + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
-                     &             ) * tmask(ji,jj,1)
-                  IF ( ln_pnd_H12 ) THEN   ! melt ponds
-                     pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps &
-                        &              + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt &
-                        &              ) * tmask(ji,jj,1)
-                  ENDIF
-              END DO
-            END DO
-            CALL lbc_lnk( pa_i(:,:,jl), 'T',  1. )
-            !
-         END DO
-
-      END DO
-      !
-   END SUBROUTINE ice_dyn_adv_umx
-
-   
-   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE adv_umx  ***
-      !! 
-      !! **  Purpose :   Compute the now trend due to total advection of 
-      !!       tracers and add it to the general trend of tracer equations
-      !!
-      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
-      !!       corrected flux (monotonic correction)
-      !!       note: - this advection scheme needs a leap-frog time scheme
-      !!
-      !! ** Action : - pt  the after advective tracer
-      !!----------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0)
-      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2)
-      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration
-      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration
-      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt             ! tracer field
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pfu, pfv       ! high order fluxes
-      !
-      INTEGER  ::   ji, jj           ! dummy loop indices  
-      REAL(wp) ::   ztra             ! local scalar
-      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt, zz1_a
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc 
-      !!----------------------------------------------------------------------
-      !
-      !  upstream (_ups) advection with initial mass fluxes
-      ! ---------------------------------------------------
-      IF( .NOT. ll_upsxy ) THEN
-
-         ! fluxes
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1
-               zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
-               zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
-            END DO
-         END DO
-
-      ELSE
-         !
-         IF    ( pamsk == 1. ) THEN ; zz1_a(:,:) = 1._wp     ;        ! 1   if advection of A
-         ELSEIF( pamsk == 0. ) THEN ; zz1_a(:,:) = z1_a(:,:) ; ENDIF  ! 1/A if advection of V
-         !
-         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
-            ! flux in x-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
-               END DO
-            END DO
-            ! first guess of tracer content from u-flux
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj)
-               END DO
-            END DO
-            CALL lbc_lnk( zpt, 'T', 1. )
-            !
-            ! flux in y-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
-               END DO
-            END DO
-            !
-         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
-            ! flux in y-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
-               END DO
-            END DO
-            ! first guess of tracer content from v-flux
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj) 
-               END DO
-            END DO
-            CALL lbc_lnk( zpt, 'T', 1. )
-            !
-            ! flux in x-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
-               END DO
-            END DO
-            !
-         ENDIF
-         
-      ENDIF
-      ! guess after content field with upstream scheme
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
-               &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj)
-            zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      CALL lbc_lnk( zt_ups, 'T', 1. )
-      
-      ! High order (_ho) fluxes 
-      ! -----------------------
-      SELECT CASE( kn_umx )
-         !
-      CASE ( 20 )                          !== centered second order ==!
-         !
-         CALL cen2( 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
-         
-      ! output high order fluxes
-      ! ------------------------
-      IF( PRESENT(pfu) ) THEN 
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1
-               pfu(ji,jj) =  zfu_ho(ji,jj)
-               pfv(ji,jj) =  zfv_ho(ji,jj)
-            END DO
-         END DO
-         !!CALL lbc_lnk( pfu, 'U',  -1. ) ! clem: not needed I think
-         !!CALL lbc_lnk( pfv, 'V',  -1. )
-      ENDIF
-      
-      ! final trend with corrected fluxes
-      ! ------------------------------------
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1 
-            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
-               &         ) * r1_e1e2t(ji,jj)  
-            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      CALL lbc_lnk( ptc, 'T',  1. )
-      !
-   END SUBROUTINE adv_umx
-
-   SUBROUTINE cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
-      &             pt_ups, pfu_ups, pfv_ups )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE macho  ***
-      !!     
-      !! **  Purpose :   compute  
-      !!
-      !! **  Method  :   ... ???
-      !!                 TIM = transient interpolation Modeling 
-      !!
-      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
-      !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
-      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
-      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
-      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes 
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      LOGICAL  ::   ll_xy = .TRUE.
-      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
-      !!----------------------------------------------------------------------
-      !
-      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --!
-         !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1
-               pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
-               pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
-            END DO
-         END DO
-         IF    ( kn_limiter == 1 ) THEN
-            CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
-         ELSEIF( kn_limiter == 2 ) THEN
-            CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
-            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
-         ELSEIF( kn_limiter == 3 ) THEN
-            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
-            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
-         ENDIF
-         !
-      ELSE                    !-- alternate directions --!
-         !
-         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
-            !
-            ! flux in x-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
-               END DO
-            END DO
-            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
-            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
-
-            ! first guess of tracer content from u-flux
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
-               END DO
-            END DO
-            CALL lbc_lnk( zzt, 'T', 1. )
-
-            ! flux in y-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
-               END DO
-            END DO
-            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
-            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
-
-         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
-            !
-            ! flux in y-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
-               END DO
-            END DO
-            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
-            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
-            !
-            ! first guess of tracer content from v-flux
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
-               END DO
-            END DO
-            CALL lbc_lnk( zzt, 'T', 1. )
-            !
-            ! flux in x-direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1
-                  pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
-               END DO
-            END DO
-            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
-            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
-
-         ENDIF
-         IF( kn_limiter == 1 )   CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
-         
-      ENDIF
-   
-   END SUBROUTINE cen2
-
-   
-   SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, &
-      &              pt_ups, pfu_ups, pfv_ups )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE macho  ***
-      !!     
-      !! **  Purpose :   compute  
-      !!
-      !! **  Method  :   ... ???
-      !!                 TIM = transient interpolation Modeling 
-      !!
-      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
-      !!----------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
-      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
-      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
-      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
-      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
-      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes 
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   ztra
-      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho
-      !!----------------------------------------------------------------------
-      !
-      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
-         !
-         !                                                        !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
-         !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
-         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
-         !                                                        !--  advective form update in zzt  --!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
-                  &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
-            END DO
-         END DO
-         CALL lbc_lnk( zzt, 'T', 1. )
-         !                                                        !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
-         !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
-         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
-         !
-      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
-         !
-         !                                                        !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
-         !                                                        !--  limiter in y --!
-         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
-         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
-         !                                                        !--  advective form update in zzt  --!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
-                  &                    - pt  (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
-            END DO
-         END DO
-         CALL lbc_lnk( zzt, 'T', 1. )
-         !                                                        !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
-         !                                                        !--  limiter in x --!
-         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
-         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
-         !
-      ENDIF
-      IF( kn_limiter == 1 ) THEN
-         IF( .NOT. ll_limiter_it2 ) THEN
-            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
-         ELSE
-            zzfu_ho(:,:) = pfu_ho(:,:)
-            zzfv_ho(:,:) = pfv_ho(:,:)
-            ! 1st iteration of nonosc (limit the flux with the upstream solution)
-            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
-            ! guess after content field with high order
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)
-                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
-               END DO
-            END DO
-            CALL lbc_lnk( zzt, 'T', 1. )
-            ! 2nd iteration of nonosc (limit the flux with the limited high order solution)
-            CALL nonosc_2d ( pdt, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
-         ENDIF
-      ENDIF
-      !
-   END SUBROUTINE macho
-
-
-   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
-      !!---------------------------------------------------------------------
-      !!                    ***  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   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
-      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity * A component
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho    ! high order flux 
-      !
-      INTEGER  ::   ji, jj             ! dummy loop indices
-      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztu1, ztu2, ztu3, ztu4
-      !!----------------------------------------------------------------------
-      !
-      !                                                     !--  Laplacian in i-direction  --!
-      DO jj = 2, jpjm1         ! First derivative (gradient)
-         DO ji = 1, fs_jpim1
-            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
-         END DO
-         !                     ! Second derivative (Laplacian)
-         DO ji = fs_2, fs_jpim1
-            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
-         END DO
-      END DO
-      CALL lbc_lnk( ztu2, 'T', 1. )
-      !
-      !                                                     !--  BiLaplacian in i-direction  --!
-      DO jj = 2, jpjm1         ! Third derivative
-         DO ji = 1, fs_jpim1
-            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
-         END DO
-         !                     ! Fourth derivative
-         DO ji = fs_2, fs_jpim1
-            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
-         END DO
-      END DO
-      CALL lbc_lnk( ztu4, 'T', 1. )
-      !
-      !
-      SELECT CASE (kn_umx )
-      !
-      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
-         !        
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
-                  &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
-            END DO
-         END DO
-         !
-      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
-         !
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
-               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
-                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
-            END DO
-         END DO
-         !  
-      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
-         !
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        &
-                  &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
-                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        &
-                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
-            END DO
-         END DO
-         !
-      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
-         !
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        &
-                  &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
-                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        &
-                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
-            END DO
-         END DO
-         !
-      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
-         !
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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)
-               zdx4 = zdx2 * zdx2
-               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       &
-                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   &
-                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       &
-                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   &
-                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       &
-                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
-            END DO
-         END DO
-         !
-      END SELECT
-      !                                                     !-- High order flux in i-direction  --!
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
-         END DO
-      END DO
-      !
-   END SUBROUTINE ultimate_x
-   
- 
-   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE ultimate_y  ***
-      !!     
-      !! **  Purpose :   compute  
-      !!
-      !! **  Method  :   ... ???
-      !!                 TIM = transient interpolation Modeling 
-      !!
-      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
-      !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
-      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity*A component
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfv_ho    ! high order flux 
-      !
-      INTEGER  ::   ji, jj       ! dummy loop indices
-      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztv1, ztv2, ztv3, ztv4
-      !!----------------------------------------------------------------------
-      !
-      !                                                     !--  Laplacian in j-direction  --!
-      DO jj = 1, jpjm1         ! First derivative (gradient)
-         DO ji = fs_2, fs_jpim1
-            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * 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) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
-         END DO
-      END DO
-      CALL lbc_lnk( ztv2, 'T', 1. )
-      !
-      !                                                     !--  BiLaplacian in j-direction  --!
-      DO jj = 1, jpjm1         ! First derivative
-         DO ji = fs_2, fs_jpim1
-            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * 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) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
-         END DO
-      END DO
-      CALL lbc_lnk( ztv4, 'T', 1. )
-      !
-      !
-      SELECT CASE (kn_umx )
-      !
-      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
-                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
-            END DO
-         END DO
-         !
-      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
-                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
-            END DO
-         END DO
-         CALL lbc_lnk( pt_v, 'V',  1. )
-         !
-      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, 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)
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       &
-                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
-                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       &
-                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
-            END DO
-         END DO
-         !
-      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, 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)
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       &
-                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
-                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       &
-                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
-            END DO
-         END DO
-         !
-      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, 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)
-               zdy4 = zdy2 * zdy2
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      &
-                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  &
-                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      &
-                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  &
-                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      &
-                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
-            END DO
-         END DO
-         !
-      END SELECT
-      !                                                     !-- High order flux in j-direction  --!
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
-         END DO
-      END DO
-      !
-   END SUBROUTINE ultimate_y
-     
-
-   SUBROUTINE nonosc_2d( pdt, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE nonosc  ***
-      !!     
-       !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
-      !!       scheme and the before field by a nonoscillatory algorithm 
-      !!
-      !! **  Method  :   ... ???
-      !!       warning : ptc and pt_low must be masked, but the boundaries
-      !!       conditions on the fluxes are not necessary zalezak (1979)
-      !!       drange (1995) multi-dimensional forward-in-time and upstream-
-      !!       in-space based differencing for fluid
-      !!----------------------------------------------------------------------
-      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_low      ! before field & upstream guess of after field
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_low, pfu_low ! upstream flux
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt   ! local scalars
-      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv, zti_low, ztj_low
-      !!----------------------------------------------------------------------
-      zbig = 1.e+40_wp
-      zsml = 1.e-15_wp
-
-      ! test on divergence
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.  
-            zdiv(ji,jj) =  - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) )  
-         END DO
-      END DO
-      CALL lbc_lnk( zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
-
-      ! antidiffusive flux : high order minus low order
-      ! --------------------------------------------------
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj)
-            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj)
-         END DO
-      END DO
-
-      ! extreme case where pfu_ho has to be zero
-      ! ----------------------------------------
-      !                                    pfu_ho
-      !                           *         --->
-      !                        |      |  *   |        | 
-      !                        |      |      |    *   |    
-      !                        |      |      |        |    *
-      !            t_low :       i-1     i       i+1       i+2   
-      IF( ll_prelimiter ) THEN
-
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1 
-               zti_low(ji,jj)= pt_low(ji+1,jj  )
-               ztj_low(ji,jj)= pt_low(ji  ,jj+1)
-            END DO
-         END DO
-         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
-
-         IF( ll_prelimiter_zalesak ) THEN
-            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) ) ) THEN
-                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj) - zti_low(ji  ,jj) ) )   pfu_ho(ji,jj) = 0.
-                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj) ) )   pfu_ho(ji,jj) = 0.
-                  ENDIF
-                  IF(    SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj+1) - pt_low (ji,jj  ) ) ) THEN
-                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) )   pfv_ho(ji,jj) = 0.
-                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj  ) - pt_low (ji,jj-1) ) )   pfv_ho(ji,jj) = 0.
-                  ENDIF
-               END DO
-            END DO
-            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
-
-         ELSEIF( ll_prelimiter_devore ) THEN
-            z1_dt = 1._wp / pdt
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) )
-                  pfu_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , &
-                     &                          zsign * ( pt_low (ji  ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji  ,jj) * z1_dt , &
-                     &                          zsign * ( zti_low(ji+1,jj) - zti_low(ji  ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) )
-                  
-                  zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) )
-                  pfv_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , &
-                     &                          zsign * ( pt_low (ji,jj  ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj  ) * z1_dt , &
-                     &                          zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) * e1e2t(ji,jj+1) * z1_dt ) )
-               END DO
-            END DO
-            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
-            
-         ENDIF
-         
-      ENDIF
-      
-      ! Search local extrema
-      ! --------------------
-      ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN
-               zbup(ji,jj) = -zbig
-               zbdo(ji,jj) =  zbig
-            ELSE               
-               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) )
-               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) )
-            ENDIF
-         END DO
-      END DO
-
-      z1_dt = 1._wp / pdt
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            !
-            !zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), 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) )
-            !
-            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)  )
-            !
-            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux
-               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )
-            zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) &
-               & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) )
-            !
-            !                                  ! up & down beta terms
-!!clem            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
-!!clem            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
-            IF( zpos >= epsi20 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
-            ELSE                      ; zbetup(ji,jj) = zbig
-            ENDIF
-            !
-            IF( zneg >= epsi20 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
-            ELSE                      ; zbetdo(ji,jj) = zbig
-            ENDIF
-            !
-         END DO
-      END DO
-      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
-
-      ! monotonic flux in the y direction
-      ! ---------------------------------
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
-            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
-            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
-            !
-            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj)
-         END DO
-      END DO
-
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
-            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
-            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
-            !
-            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj)
-         END DO
-      END DO
-      !
-   END SUBROUTINE nonosc_2d
-
-   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE limiter_x  ***
-      !!     
-      !! **  Purpose :   compute flux limiter 
-      !!----------------------------------------------------------------------
-      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfu_ho       ! high order flux
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux
-      !
-      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp), DIMENSION (jpi,jpj) ::   zslpx       ! tracer slopes 
-      !!----------------------------------------------------------------------
-      !
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1)
-         END DO
-      END DO
-      CALL lbc_lnk( zslpx, 'U', -1.)   ! lateral boundary cond.
-      
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
-
-            Rjm = zslpx(ji-1,jj)
-            Rj  = zslpx(ji  ,jj)
-            Rjp = zslpx(ji+1,jj)
-
-            IF( PRESENT(pfu_ups) ) THEN
-
-               IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
-               ELSE                        ;   Rr = Rjp
-               ENDIF
-
-               zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj)     
-               IF( Rj > 0. ) THEN
-                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)),  &
-                     &        MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
-               ELSE
-                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)),  &
-                     &        MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
-               ENDIF
-               pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter
-               
-            ELSE
-               IF( Rj /= 0. ) THEN
-                  IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
-                  ELSE                        ;   Cr = Rjp / Rj
-                  ENDIF
-               ELSE
-                  Cr = 0.
-                  !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
-                  !ELSE                        ;   Cr = Rjp * 1.e20
-                  !ENDIF
-               ENDIF
-               
-               ! -- superbee --
-               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
-               ! -- van albada 2 --
-               !!zpsi = 2.*Cr / (Cr*Cr+1.)
-
-               ! -- sweby (with beta=1) --
-               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
-               ! -- van Leer --
-               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
-               ! -- ospre --
-               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
-               ! -- koren --
-               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
-               ! -- charm --
-               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
-               !ELSE                 ;   zpsi = 0.
-               !ENDIF
-               ! -- van albada 1 --
-               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
-               ! -- smart --
-               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
-               ! -- umist --
-               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
-
-               ! high order flux corrected by the limiter
-               pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
-
-            ENDIF
-         END DO
-      END DO
-      CALL lbc_lnk( pfu_ho, 'U', -1.)   ! lateral boundary cond.
-      !
-   END SUBROUTINE limiter_x
-
-   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE limiter_y  ***
-      !!     
-      !! **  Purpose :   compute flux limiter 
-      !!----------------------------------------------------------------------
-      REAL(wp)                     , INTENT(in   )           ::   pdt          ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   )           ::   pt           ! ice tracer
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout)           ::   pfv_ho       ! high order flux
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux
-      !
-      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp), DIMENSION (jpi,jpj) ::   zslpy       ! tracer slopes 
-      !!----------------------------------------------------------------------
-      !
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1)
-         END DO
-      END DO
-      CALL lbc_lnk( zslpy, 'V', -1.)   ! lateral boundary cond.
-      
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
-
-            Rjm = zslpy(ji,jj-1)
-            Rj  = zslpy(ji,jj  )
-            Rjp = zslpy(ji,jj+1)
-
-            IF( PRESENT(pfv_ups) ) THEN
-
-               IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
-               ELSE                        ;   Rr = Rjp
-               ENDIF
-
-               zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj)     
-               IF( Rj > 0. ) THEN
-                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)),  &
-                     &        MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
-               ELSE
-                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)),  &
-                     &        MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
-               ENDIF
-               pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter
-
-            ELSE
-
-               IF( Rj /= 0. ) THEN
-                  IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
-                  ELSE                        ;   Cr = Rjp / Rj
-                  ENDIF
-               ELSE
-                  Cr = 0.
-                  !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
-                  !ELSE                        ;   Cr = Rjp * 1.e20
-                  !ENDIF
-               ENDIF
-
-               ! -- superbee --
-               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
-               ! -- van albada 2 --
-               !!zpsi = 2.*Cr / (Cr*Cr+1.)
-
-               ! -- sweby (with beta=1) --
-               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
-               ! -- van Leer --
-               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
-               ! -- ospre --
-               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
-               ! -- koren --
-               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
-               ! -- charm --
-               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
-               !ELSE                 ;   zpsi = 0.
-               !ENDIF
-               ! -- van albada 1 --
-               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
-               ! -- smart --
-               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
-               ! -- umist --
-               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
-
-               ! high order flux corrected by the limiter
-               pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
-
-            ENDIF
-         END DO
-      END DO
-      CALL lbc_lnk( pfv_ho, 'V', -1.)   ! lateral boundary cond.
-      !
-   END SUBROUTINE limiter_y
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option           Dummy module         NO SI3 sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE icedyn_adv_umx
Index: NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEDYN/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEDYN/EXPREF/namelist_ice_cfg	(revision 10397)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEDYN/EXPREF/namelist_ice_cfg	(revision 10399)
@@ -33,7 +33,8 @@
 &namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
-   ln_dynFULL       = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   ln_dynALL        = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
    ln_dynRHGADV     = .true.          !  dyn.: no ridge/raft & no corrections  (rheology + advection)
-   ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection)
+   ln_dynADV1D      = .false.         !  dyn.: only advection 1D                  (Schar & Smolarkiewicz 1996 test case)
+   ln_dynADV2D      = .false.         !  dyn.: only advection 2D w prescribed vel.(rn_uvice + advection)
       rn_uice       =   0.5           !        prescribed ice u-velocity
       rn_vice       =   0.            !        prescribed ice v-velocity
