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 10315)
+++ NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90	(revision 10331)
@@ -37,4 +37,7 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a, z1_a_ups, zua_ups, zva_ups
    
+   ! rachid trick (only for nonosc)
+   LOGICAL :: ll_rachid = .TRUE.
+
    ! alternate directions for upstream
    ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false"
@@ -44,7 +47,6 @@
    ! 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
+   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11
 
    ! iterate on the limiter (only nonosc)
@@ -154,10 +156,15 @@
 
             ! 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) )
+!!            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
@@ -193,4 +200,5 @@
                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
@@ -256,4 +264,5 @@
       !  upstream (_ups) advection with initial mass fluxes
       ! ---------------------------------------------------
+      IF( ll_rachid )    zfu_ups=0.; zfv_ups=0.
       IF( .NOT. ll_upsxy ) THEN
 
@@ -320,4 +329,23 @@
          
       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
+      ENDIF
+
       ! guess after content field with upstream scheme
       DO jj = 2, jpjm1
@@ -325,5 +353,9 @@
             ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
                &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj)
-            zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+            IF( ll_rachid ) THEN
+               zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+            ELSE
+               zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+            ENDIF
          END DO
       END DO
@@ -336,5 +368,5 @@
       CASE ( 20 )                          !== centered second order ==!
          !
-         CALL cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
+         CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
             &       zt_ups, zfu_ups, zfv_ups )
          !
@@ -377,5 +409,5 @@
    END SUBROUTINE adv_umx
 
-   SUBROUTINE cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
+   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
       &             pt_ups, pfu_ups, pfv_ups )
       !!---------------------------------------------------------------------
@@ -389,4 +421,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
       INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
       INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
@@ -399,5 +432,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes 
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
@@ -415,5 +448,5 @@
          END DO
          IF    ( kn_limiter == 1 ) THEN
-            CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+            CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
          ELSEIF( kn_limiter == 2 ) THEN
             CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
@@ -483,5 +516,5 @@
 
          ENDIF
-         IF( kn_limiter == 1 )   CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
          
       ENDIF
@@ -516,5 +549,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes 
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pfu_ups, pfv_ups ! upstream fluxes 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes 
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
@@ -570,10 +603,10 @@
       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 )
+            CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, 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 )
+            CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
             ! guess after content field with high order
             DO jj = 2, jpjm1
@@ -585,5 +618,5 @@
             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 )
+            CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
          ENDIF
       ENDIF
@@ -842,5 +875,5 @@
      
 
-   SUBROUTINE nonosc_2d( pdt, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
+   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -855,7 +888,12 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
+      REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
       REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   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   ) ::   pfv_low, pfu_low ! upstream flux
+      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
       !
@@ -868,4 +906,18 @@
       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         
+      ENDIF
+      ENDIF
+      
       ! antidiffusive flux : high order minus low order
       ! --------------------------------------------------
@@ -885,6 +937,6 @@
       !                        |      |      |        |    *
       !            t_low :       i-1     i       i+1       i+2   
-      IF( ll_prelimiter ) THEN
-
+      IF( ll_prelimiter_zalesak ) THEN
+         
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1 
@@ -895,62 +947,68 @@
          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.
+
+         !! 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
+         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. )
+
+         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
          
-      ENDIF
       
       ! Search local extrema
@@ -982,6 +1040,6 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             !
-            !zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1) )  ! search max/min in neighbourhood
-            !zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) )
+!            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1) )  ! search max/min in neighbourhood
+!            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
@@ -1014,4 +1072,26 @@
       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
       ! ---------------------------------
