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 10312)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90	(revision 10315)
@@ -35,4 +35,22 @@
    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
 
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a, z1_a_ups, zua_ups, zva_ups
+   
+   ! 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.
+   
+   ! prelimiter
+   ! clem: use it to avoid overshoot in H
+   LOGICAL :: ll_prelimiter = .TRUE.
+   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
+
+   ! iterate on the limiter (only nonosc)
+   ! clem: useless
+   LOGICAL :: ll_limiter_it2 = .FALSE.
+   
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -74,8 +92,8 @@
       REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
       REAL(wp) ::   zcfl , zdt
-      REAL(wp) ::   zeps = 0.1_wp           ! shift in concentration to avoid division by 0
+      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) ::   z1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip 
+      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 
       !!----------------------------------------------------------------------
       !
@@ -112,4 +130,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))
       !---------------!
       !== advection ==!
@@ -117,6 +139,6 @@
       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 
+         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 
 
          DO jl = 1, jpl
@@ -125,40 +147,52 @@
             !
             WHERE( pa_i(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_i(:,:,jl)
-            ELSEWHERE                        ;   z1_a(:,:) = 0.   !!;   pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0. 
+            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_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 = 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
+            WHERE( za_ups(:,:) > epsi20 )   ;   z1_a_ups(:,:) = 1._wp / za_ups(:,:)
+            ELSEWHERE                       ;   z1_a_ups(:,:) = 0.
+            END WHERE
+
+            !!zua_ho = zudy ; zva_ho = zvdx
+            !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zua_ho , zva_ho , 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
+            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, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) )              ! Snw volume
+            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, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) )              ! Salt content
+            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, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl) )              ! Age content
+            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
             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
+               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
             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
+               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)
                !
-               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_ap(:,:) = 1._wp / pa_ip(:,:,jl)
-               ELSEWHERE                         ;   z1_ap(:,:) = 0.
+               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 = 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
+
                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
+               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
             !
@@ -177,5 +211,5 @@
               END DO
             END DO
-            CALL lbc_lnk( pa_i(:,:,jl), 'T',  1. )
+            CALL lbc_lnk_multi( pa_i(:,:,jl), 'T',  1., pa_ip(:,:,jl), 'T',  1. )
             !
          END DO
@@ -186,5 +220,5 @@
 
    
-   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv )
+   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 )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -209,9 +243,10 @@
       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
+      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
       !
       INTEGER  ::   ji, jj           ! dummy loop indices  
       REAL(wp) ::   ztra             ! local scalar
-!!clem      REAL(wp) ::   zeps = 1.e-02        ! local scalar
       INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt
@@ -219,20 +254,70 @@
       !!----------------------------------------------------------------------
       !
-      ! add a constant value to avoid problems with zeros
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1)
-         END DO
-      END DO
-
       !  upstream (_ups) advection with initial mass fluxes
       ! ---------------------------------------------------
-      ! fluxes
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1
-            zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
-            zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
-         END DO
-      END DO
+      IF( .NOT. ll_upsxy ) THEN
+
+         ! fluxes
+         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)
+            END DO
+         END DO
+
+      ELSE
+         ! 1 if advection of A
+         ! z1_a_ups already defined IF advection of other tracers
+         IF( pamsk == 1. )   z1_a_ups(:,:) = 1._wp
+         !
+         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+            ! flux in x-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  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)
+               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)
+               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( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pva_ups(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( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(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) * z1_a_ups(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( pua_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pua_ups(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
@@ -244,5 +329,5 @@
       END DO
       CALL lbc_lnk( zt_ups, 'T', 1. )
-      
+
       ! High order (_ho) fluxes 
       ! -----------------------
@@ -251,25 +336,30 @@
       CASE ( 20 )                          !== centered second order ==!
          !
-         CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
+         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, zpt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
+         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 
+      ! 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
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1
-               pfu(ji,jj) =  zfu_ho(ji,jj) !!clem - zeps * puc(ji,jj)
-               pfv(ji,jj) =  zfv_ho(ji,jj) !!clem - zeps * pvc(ji,jj)
-            END DO
-         END DO
-         !!CALL lbc_lnk( pfu, 'U',  -1. ) ! clem: not needed I think
-         !!CALL lbc_lnk( pfv, 'V',  -1. )
+               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)
+            END DO
+         END DO
+         !!CALL lbc_lnk( pua_ho, 'U',  -1. ) ! clem: not needed I think
+         !!CALL lbc_lnk( pva_ho, 'V',  -1. )
       ENDIF
       
@@ -278,6 +368,5 @@
       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)
-!!clem               &           + ( puc   (ji,jj) - puc   (ji-1,jj) + pvc   (ji,jj) - pvc   (ji,jj-1) ) * zeps  &  ! epsi * Div(ua)  or epsi * Div(u) 
+            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)
@@ -430,5 +519,6 @@
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
+      REAL(wp) ::   ztra
+      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho
       !!----------------------------------------------------------------------
       !
@@ -478,5 +568,24 @@
          !
       ENDIF
-      IF( kn_limiter == 1 )   CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+      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
@@ -733,13 +842,13 @@
      
 
-   SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+   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 
+       !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
       !!       scheme and the before field by a nonoscillatory algorithm 
       !!
       !! **  Method  :   ... ???
-      !!       warning : ptc and pt_ups must be masked, but the boundaries
+      !!       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-
@@ -747,23 +856,15 @@
       !!----------------------------------------------------------------------
       REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_ups      ! before field & upstream guess of after field
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(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         !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, 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)
+      zsml = epsi20
 
       ! antidiffusive flux : high order minus low order
@@ -771,20 +872,106 @@
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector opt.
-            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_ups(ji,jj)
-            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_ups(ji,jj)
-         END DO
-      END DO
-
+            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
+
+            !! this does not work
+!!            DO jj = 2, jpjm1
+!!               DO ji = fs_2, fs_jpim1
+!!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj  ) - pt_low (ji  ,jj) ) .AND.     &
+!!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj+1) - pt_low (ji  ,jj) )           &
+!!               &    ) THEN
+!!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji  ,jj) ) .AND.   &
+!!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji  ,jj) )         &
+!!               &       ) THEN
+!!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0.
+!!                     ENDIF
+!!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj  ) ) .AND.  &
+!!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji  ,jj-1) )        &
+!!               &       )  THEN
+!!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0.
+!!                     ENDIF
+!!                  ENDIF
+!!                END DO
+!!            END DO            
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND.  &
+                     & pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN
+                     !
+                     IF(  pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND.  &
+                        & pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0.
+                     !
+                     IF(  pfu_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji-1,jj) ) <= 0 .AND.  &
+                        & pfv_ho(ji,jj) * ( pt_low(ji  ,jj) - pt_low(ji,jj-1) ) <= 0)  pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0.
+                     !
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+
+         ELSEIF( ll_prelimiter_devore ) THEN
+            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_ups with large negative/positive value (-/+zbig) outside ice cover
+      ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover
       DO jj = 1, jpj
-         DO ji = fs_2, fs_jpim1
-            IF( ptc(ji,jj) == 0._wp .AND. pt_ups(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN
+         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
                zbup(ji,jj) = -zbig
                zbdo(ji,jj) =  zbig
-            ELSE               
-               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_ups(ji,jj) )
-               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_ups(ji,jj) )
+!!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) = 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)
+            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
@@ -795,26 +982,31 @@
          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) )
-            !
-            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &
-               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )  ! positive/negative part of the flux
+            !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_ups(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
-!!clem            zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
-            IF( zpos >= epsi20 ) THEN
-               zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
-            ELSE
-               zbetup(ji,jj) = zbig
+!            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 >= epsi20 ) THEN
-               zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
-            ELSE
-               zbetdo(ji,jj) = zbig
+            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.
+            !
             !
          END DO
@@ -830,5 +1022,5 @@
             zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
             !
-            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj)
+            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj)
          END DO
       END DO
@@ -840,5 +1032,5 @@
             zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
             !
-            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj)
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj)
          END DO
       END DO
