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 10331)
+++ 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
